{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from genomics.popgen.admix import cluster, plot\n",
    "\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "k_range = range(2, 10)  # 2..9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The next cell is very slow. Example outputs are provided (so you can avoid running it)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#for k in k_range:\n",
    "#    os.system('admixture --cv=10 hapmap10_auto_noofs_ld.bed %d > admix.%d' % (k, k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Individual order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = open('hapmap10_auto_noofs_ld.fam')\n",
    "ind_order = []\n",
    "for l in f:\n",
    "    toks = l.rstrip().replace(' ', '\\t').split('\\t')\n",
    "    fam_id = toks[0]\n",
    "    ind_id = toks[1]\n",
    "    ind_order.append((fam_id, ind_id))\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CV-plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'K')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7IAAAImCAYAAABuALYVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3WmQXWd+Hvbn7W7sG4mtbw83cCfQjVkhamY0O4dAN+yM5FjlSJmSbKdilVxWVLFUthQviiIldqLEkuOKbJfssmNlIsuy4tgTGQs5C2ffMAsH3QA5BDnkECS6sRLEvnS/+dCXNAYCiQbZwOnb/ftV3QLvuadPPwf4MPXM+d/3LbXWAAAAQKfoajoAAAAAXAtFFgAAgI6iyAIAANBRFFkAAAA6iiILAABAR1FkAQAA6CiKLADcIKWUv1RK+eIl70+WUu6ayrlv4HdtL6X8xTf68wAwkymyAMxYpZT/spSyq134DrTL2fsayrKwlPJSKeUjV/jsd0opf3yt16y1Lq21PjMN2X69lPKJy649VGv9V2/22gAwEymyAMxIpZRfSvIPk/y9JL1Jbk/yj5P8+Guc33M989Razyb5N0l+9rLf253kp5MojZe50r/Jtf47Xe9/VwA6kyILwIxTSlmR5DeS/LVa67+rtZ6qtV6otf5/tda/0T7n10spf1xK+UQp5eUkf6mUsqCU8g9LKS+2X/+wlLKgff7qUsqftJ+qHi2lfKGU0tX+7FdKKS+UUk6UUp4spTz0GtH+VZI/X0pZfMmxLZn839Pt7Wv9ainl6fa19pRS/tzr3GctpdzT/u9VpZRPllJeLqV8Pcndl537v5dSnm9//s1SyvvbxweT/K0k/0X7yfXj7eOPlVL+6/Z/d5VS/k4p5blSysFSyu+3/45TSlnXzvEXSyk/KKUcLqX87dfJvKCU8r+1zx0rpfzTUsqi9mcfKqXsb/99jib5l1c61j73r5RS9rX/LT5ZSnnLZX8vf62U8lSSp14rCwBzlyILwEz0niQLk/y/Vznvx5P8cZKbkvzfSf52kncneXuStyV5MMnfaZ/7y0n2J1mTySe8fytJLaXcn+QXkvxIrXVZJovps1f6ZbXWLyc5kOQ/v+TwzyT5g1rrxfb7p5O8P8mKJP9Dkk+UUvqmcM+/m+Rskr4k/1X7dalvtO9rZZI/SPJvSykLa607MvnU+t+0R5XfdoVr/6X268NJ7kqyNMn/cdk570tyf5KHkvxaKWX9a+T8X5Lc185yT5JbkvzaJZ+32hnvSPJzVzrWHs/++0n+Qvt+n0vyh5f9np9I8qNJNrxGDgDmMEUWgJloVZLDl5TD1/KVWuu/r7VO1FrPJPl4kt+otR6stR7KZJH8mfa5FzJZmu5oP939Qq21JhlPsiDJhlLKvFrrs7XWp1/nd/5+2uPFpZTlmSzTr44V11r/ba31xXamf5PJJ4oPvt5NtMeT/3ySX2s/fR7OZaPKtdZP1FqP1Fov1lr/QTvz/Vf5+3nFx5P8dq31mVrryST/XZKfumxs93+otZ6ptT6e5PFM/h8Bl+csSf5Kkr9eaz1aaz2RyRL9U5ecNpHkv6+1nmv/m1zp2MeT/Ita67dqrefaed5TSll3yXX+fvt3nAkAXEaRBWAmOpJk9RS+H/n8Ze/fksmne694rn0sSf7XJPuSPFJKeaaU8qtJUmvdl+S/TfLrSQ6WUv7wlTHX9qjuK6/b29f5/SQfLqXckuQnk+yrtX77lV9YSvnZUsp32iPMLyUZSLL6KvexJknPZfdz6X2klPLLpZS9pZTj7euumMJ1X3Glv5eeTD6ZfsXoJf99OpNPba+Uc3GSb15yfzvax19xqP194rzOsR/K0y7XRzL5dPcVl//bAsCrFFkAZqKvZHLM9ieucl697P2LmRxffcXt7WOptZ6otf5yrfWuJP9Zkl965buwtdY/qLW+r/2zNZPjs6+sKvzK6wftYz9I8oVMPlX8mUwW2yRJKeWOJP8sk6PKq2qtNyUZTlKuch+HklxMcttl2V+57vuT/EomR3Fvbl/3+CXXvfzv4XJX+nu5mGTsKj93ucNJziTpr7Xe1H6tqLVeWnqvlOV1/51KKUsy+RT+hatcBwCSKLIAzEC11uOZ/N7l75ZSfqKUsriUMq+UMlRK+a3X+dF/neTvlFLWlFJWt6/xiSQppfzZUso97fHYlzM5UjxeSrm/lPKR9qJQZzNZ1MavEvFfZbKs/lgmv5v7iiWZLGCH2r/zL2fyiezV7nc8yb9L8uvte92Q5NI9YJdlsngeStJTSvm1JMsv+XwsybpXFq+6gn+d5K+XUu4spSzNf/pO7dVGty/POZHJov47pZS17Xu8pZSy5Vquk8nv+P7lUsrb23/vfy/J12qtz17jdQCYoxRZAGakWutvJ/mlTC7WdCiTo6a/kOTfv86P/Y9JdiX5bpLdSb7VPpYk9yb5VJKTmXzi+49rrY9l8rum/3MmnzaOJlmbyYWgXs8fJ7k5yadrrQcuybwnyT9oX38sycYkX5rK/bbvbWk7w/+Z9uq+bTszuSry9zI5kns2Pzx6+2/bfx4ppXzrCtf+F0n+rySfT/L99s//N1PMdblfyeSI9lfL5GrRn8rUv6ubJKm1fjrJ303y/2Ry8ay788PfswWA11Um17kAAACAzuCJLAAAAB1FkQUAAKCjKLIAAAB0FEUWAACAjqLIAgAA0FF6mg5wLVavXl3XrVvXdAwAAACug29+85uHa61rrnZeRxXZdevWZdeuXU3HAAAA4DoopTw3lfOMFgMAANBRFFkAAAA6iiILAABAR1FkAQAA6CiKLAAAAB1FkQUAAKCjKLIAAAB0FEUWAACAjqLIAgAA0FEUWQAAADqKIgsAAEBHUWQBAADoKIosAAAAHUWRBQAAoKMosgAAAHQURRYAAICOosgCAADQURRZAAAAOooiO43OnB9vOgIAAMCsp8hOk9/7/NN599//dM5dVGYBAACuJ0V2mtzXuyzHz1zIF5863HQUAACAWU2RnSbvvXt1li3sybbdo01HAQAAmNUU2Wkyv6crD2/ozaf2juXC+ETTcQAAAGYtRXYaDQ305fiZC/nK00eajgIAADBrKbLT6P33rs6S+d3ZPnyg6SgAAACzliI7jRbO685H1vfmkZGxXDReDAAAcF0ostNsaKCVI6fO5+vPHm06CgAAwKykyE6zD92/JgvndWXHsNWLAQAArgdFdpotnt+TD923NjuGRzMxUZuOAwAAMOsostfB0MZWDp44l2/94FjTUQAAAGYdRfY6+MgDazO/uyvbjRcDAABMO0X2Oli2cF7ef+/q7BgeTa3GiwEAAKaTInudDG3sywsvncl39x9vOgoAAMCsosheJw+v701PVzFeDAAAMM0U2etkxeJ5ec/dq7J9+IDxYgAAgGmkyF5HWzf25bkjp7P3wImmowAAAMwaiux1tHlDb7pKsmP4QNNRAAAAZg1F9jpatXRBHrxzZbb5niwAAMC0UWSvs60b+7Lv4Mk8NWa8GAAAYDoostfZlv5Wkli9GAAAYJoostdZ7/KF2XTHzYosAADANFFkb4DBgVb2Hng5zx4+1XQUAACAjqfI3gCDA8aLAQAAposiewPcevPivO3WFbbhAQAAmAaK7A0yONCXx/cfz/5jp5uOAgAA0NEU2RtkqD1evMN4MQAAwJuiyN4g61Yvyfq+5YosAADAm6TI3kBDA63seu5Yxl4+23QUAACAjqXI3kBbN06OF+8c8VQWAADgjVJkb6B71i7LPWuXZttuqxcDAAC8UYrsDbZ1oJWvf/9ojpw813QUAACAjqTI3mCDA32ZqMkje8aajgIAANCRFNkbbH3fstyxarHxYgAAgDdIkb3BSikZGujLV54+kpdOn286DgAAQMdRZBswNNDKxYmaR40XAwAAXDNFtgFvvXVFbrlpUXYM24YHAADgWimyDSilZHCglS88dTgnzl5oOg4AAEBHUWQbMjTQyvnxiXzmiYNNRwEAAOgoimxD3nn7zVm7bEG27zZeDAAAcC0U2YZ0dU2OFz/2vYM5ff5i03EAAAA6hiLboMGBVs5emMhjTx5qOgoAAEDHmFKRLaUMllKeLKXsK6X86uuc95OllFpK2XTJsbeWUr5SShkppewupSxsH3+sfc3vtF9r3/ztdJYH163MqiXzs93qxQAAAFPWc7UTSindSX43ycNJ9if5Rinlk7XWPZedtyzJLyb52iXHepJ8IsnP1FofL6WsSnLpMr0fr7XuevO30Zl6uruyub83n/zOizl7YTwL53U3HQkAAGDGm8oT2QeT7Ku1PlNrPZ/kD5P8+BXO+80kv5Xk7CXHNif5bq318SSptR6ptY6/ycyzytBAX06dH88XnjrcdBQAAICOMJUie0uS5y95v7997FWllHckua3W+ieX/ex9SWopZWcp5VullL952ef/sj1W/HdLKeVaw88G77l7VVYsmpftwweajgIAANARplJkr1Qw66sfltKV5HeS/PIVzutJ8r4kH2//+edKKQ+1P/t4rXVjkve3Xz9zxV9eys+VUnaVUnYdOjT7FkWa192Vj67vzaN7xnL+4kTTcQAAAGa8qRTZ/Uluu+T9rUlevOT9siQDSR4rpTyb5N1JPtle8Gl/ks/VWg/XWk8n2ZbknUlSa32h/eeJJH+QyRHmP6XW+nu11k211k1r1qy5lnvrGFs3tnLi7MV8+WnjxQAAAFczlSL7jST3llLuLKXMT/JTST75yoe11uO11tW11nW11nVJvprkY+1FnHYmeWspZXF74acPJtlTSukppaxOklLKvCR/NsnwtN5ZB3nfvauzdEFPtu+2ejEAAMDVXLXI1lovJvmFTJbSvUn+qNY6Ukr5jVLKx67ys8eS/HYmy/B3knyr1vofkyxIsrOU8t328ReS/LM3dScdbEFPdz7ywNo8smc0F8eNFwMAALyeq26/kyS11m2ZHAu+9Nivvca5H7rs/ScyuQXPpcdOJXnXtQSd7bZubOWTj7+Yr3//aN57z+qm4wAAAMxYUxkt5gb44H1rs2hed7ZZvRgAAOB1KbIzxKL53fnwA2uyc2Qs4xP16j8AAAAwRymyM8jgQF8OnTiXbz53rOkoAAAAM5YiO4N85IG1md/Tle3GiwEAAF6TIjuDLF3Qkw/cuyY7hkczYbwYAADgihTZGWZooJUDx8/m8f0vNR0FAABgRlJkZ5iPru/NvO6SHcOjTUcBAACYkRTZGWbF4nl5792rs314NLUaLwYAALicIjsDDQ208oOjpzPy4stNRwEAAJhxFNkZaHN/K91dxosBAACuRJGdgVYumZ8fvXNltg0fMF4MAABwGUV2hhra2JdnDp3KUwdPNh0FAABgRlFkZ6gt/b0pJdm+23gxAADApRTZGWrtsoXZdMfN2T58oOkoAAAAM4oiO4MNDfTlidETeeaQ8WIAAIBXKLIz2OBAK0my3erFAAAAr1JkZ7C33LQob7vtJtvwAAAAXEKRneG2DrSy+4Xjef7o6aajAAAAzAiK7Aw3NNCXJJ7KAgAAtCmyM9ztqxan/y3LrV4MAADQpsh2gKGBVr71g5dy4PiZpqMAAAA0TpHtAIPt8eKdxosBAAAU2U5wz9qlua93qW14AAAAosh2jMGBvnz92aM5dOJc01EAAAAapch2iK0bW6k1eWSPp7IAAMDcpsh2iPt7l+XO1UtswwMAAMx5imyHKKVkcKCVLz99JMdOnW86DgAAQGMU2Q6ydaAv4xM1j+4dazoKAABAYxTZDjJwy/LcevOibN99oOkoAAAAjVFkO0gpJUMDrXxx3+G8fPZC03EAAAAaoch2mMGBvlwYr/nM3oNNRwEAAGiEItth3nHbTeldviDbjBcDAABzlCLbYbq6SoYG+vK57x3KqXMXm44DAABwwymyHWhwoJVzFyfy2SeNFwMAAHOPItuBfmTdyqxeOj/bh0ebjgIAAHDDKbIdqLurZHN/K5994mDOXhhvOg4AAMANpch2qKGBVk6fH8/nvneo6SgAAAA3lCLbod5916rctHhedhgvBgAA5hhFtkPN6+7Kw+t786k9Yzl30XgxAAAwdyiyHWxoYysnzl3Ml/cdaToKAADADaPIdrAfu2d1li3oyfbhA01HAQAAuGEU2Q62oKc7D61fm0f2jOXC+ETTcQAAAG4IRbbDDW3sy0unL+RrzxxtOgoAAMANoch2uA/etyaL53cbLwYAAOYMRbbDLZzXnQ/fvzY7R0YzPlGbjgMAAHDdKbKzwNDGVg6fPJ9dzxovBgAAZj9Fdhb48P1rs6CnK9uHR5uOAgAAcN0psrPAkgU9+cB9a7JjeDQTxosBAIBZTpGdJbZubGX05bP5zv6Xmo4CAABwXSmys8RHHujNvO6S7butXgwAAMxuiuwssWLRvLzvntXZPjyaWo0XAwAAs5ciO4sMDfRl/7EzGX7h5aajAAAAXDeK7Czy8IbedHeVbB82XgwAAMxeiuwscvOS+XnPXauMFwMAALOaIjvLDA608v3Dp/Lk2ImmowAAAFwXiuwss6W/lVKS7btHm44CAABwXSiys8yaZQvyI+tW+p4sAAAwaymys9DQQCvfGzuZpw+dbDoKAADAtFNkZ6HBgVaSZMew8WIAAGD2UWRnob4Vi/KO22/Ktt3GiwEAgNlHkZ2ltg70ZeTFl/ODI6ebjgIAADCtFNlZ6tXx4hFPZQEAgNlFkZ2lblu5OAO3LM822/AAAACzjCI7iw0N9OU7z7+UF18603QUAACAaaPIzmJDVi8GAABmIUV2FrtrzdLc37tMkQUAAGYVRXaWG9rYyjeeO5qDJ842HQUAAGBaKLKz3NBAX2pNdo6MNR0FAABgWiiys9x9vUtz15ol2TFsGx4AAGB2UGRnuVJKhgZa+eozR3P01Pmm4wAAALxpiuwcMDTQl/GJmkf3WPQJAADofIrsHND/luW5beWibLd6MQAAMAsosnPA5HhxX76073COn7nQdBwAAIA3RZGdI4YGWrkwXvPpvVYvBgAAOpsiO0e87dab0rdiYbbtNl4MAAB0NkV2jujqKtnS38rnnzqUk+cuNh0HAADgDVNk55CtG/ty/uJEPvvEwaajAAAAvGGK7BzyrjtuzuqlC7J9+EDTUQAAAN4wRXYO6e4q2dLfm88+cShnzo83HQcAAOANUWTnmK0b+3Lmwng+971DTUcBAAB4QxTZOeZH71yZmxfPM14MAAB0LEV2junp7srmDa18eu/BnLtovBgAAOg8iuwcNLixlZPnLuaLTx1uOgoAAMA1U2TnoB+7e3WWLezJ9uHRpqMAAABcsykV2VLKYCnlyVLKvlLKr77OeT9ZSqmllE2XHHtrKeUrpZSRUsruUsrC9vF3td/vK6X8o1JKefO3w1TM7+nKw+t78+iesVwYn2g6DgAAwDW5apEtpXQn+d0kQ0k2JPnpUsqGK5y3LMkvJvnaJcd6knwiyc/XWvuTfCjJhfbH/yTJzyW5t/0afDM3wrUZHGjl+JkL+crTR5qOAgAAcE2m8kT2wST7aq3P1FrPJ/nDJD9+hfN+M8lvJTl7ybHNSb5ba308SWqtR2qt46WUviTLa61fqbXWJL+f5CfezI1wbT5w35osmd9tvBgAAOg4UymytyR5/pL3+9vHXlVKeUeS22qtf3LZz96XpJZSdpZSvlVK+ZuXXHP/613zkmv/XCllVyll16FD9j6dLgvndefDD6zNIyOjGZ+oTccBAACYsqkU2St9d/XV5lNK6UryO0l++Qrn9SR5X5KPt//8c6WUh652zR86WOvv1Vo31Vo3rVmzZgpxmaqhgb4cOXU+X//+0aajAAAATNlUiuz+JLdd8v7WJC9e8n5ZkoEkj5VSnk3y7iSfbC/4tD/J52qth2utp5NsS/LO9vFbX+ea3AAfun9NFs7ryo7hA01HAQAAmLKpFNlvJLm3lHJnKWV+kp9K8slXPqy1Hq+1rq61rqu1rkvy1SQfq7XuSrIzyVtLKYvbCz99MMmeWuuBJCdKKe9ur1b8s0n+w/TeGlezZEFPPnjfmmwfHs2E8WIAAKBDXLXI1lovJvmFTJbSvUn+qNY6Ukr5jVLKx67ys8eS/HYmy/B3knyr1vof2x//1ST/PMm+JE8n2f6G74I3bOvGvhw8cS7ffv5Y01EAAACmpGcqJ9Vat2VyLPjSY7/2Gud+6LL3n8jkFjyXn7crkyPJNOgjD6zN/O6ubNs9mnfdsbLpOAAAAFc1ldFiZrFlC+flffeuzo7h0UzuhAQAADCzKbJkaKCVF146k90vHG86CgAAwFUpsuThDb3p6SrZtnu06SgAAABXpciSmxbPz3vuXpUdwweMFwMAADOeIkuSZGigL88eOZ0nRk80HQUAAOB1KbIkSTb396arJNt3H2g6CgAAwOtSZEmSrF66IA/euTLbh31PFgAAmNkUWV41NNCXpw6ezL6DxosBAICZS5HlVVv6W0mS7VYvBgAAZjBFlle1VizMu+642XgxAAAwoymy/JChgVb2HHg5zx051XQUAACAK1Jk+SGDA+3xYk9lAQCAGUqR5YfcevPivPXWFbbhAQAAZixFlj9lcKCVx/cfzwsvnWk6CgAAwJ+iyPKnDA30JUl2GC8GAABmIEWWP+XO1UvyQGuZ8WIAAGBGUmS5oq0b+/LNHxzL2Mtnm44CAADwQxRZrmhooJVak50jxosBAICZRZHliu7tXZa71yzJ9t2KLAAAMLMosrymrRv78rXvH8mRk+eajgIAAPAqRZbXNDjQykRNHtkz1nQUAACAVymyvKYNfctz+8rF2W4bHgAAYAZRZHlNpZQMbWzly/sO5/jpC03HAQAASKLIchVDA325OFHz6F7jxQAAwMygyPK63nbrirxlxcLsGD7QdBQAAIAkiixXUUrJ4EBfPv+9wzlx1ngxAADQPEWWqxra2Mr58Yl85omDTUcBAABQZLm6d91+c9YuW5AdVi8GAABmAEWWq+rqKtnS38pnnzyY0+cvNh0HAACY4xRZpmRoYytnL0zkc08eajoKAAAwxymyTMmD61Zm5ZL52Wa8GAAAaJgiy5T0dHdl84befGbvWM5eGG86DgAAMIcpskzZ0Ma+nDo/ni8+dbjpKAAAwBymyDJl77lrVZYv7Mm24QNNRwEAAOYwRZYpm9/TlY9u6M2n9ozl/MWJpuMAAABzlCLLNdk60JeXz17Ml582XgwAADRDkeWavO/e1Vkyvzs7rF4MAAA0RJHlmiyc152H1vfmkT1juThuvBgAALjxFFmu2dBAK0dPnc/Xv3+06SgAAMAcpMhyzT54/5osnNeV7caLAQCABiiyXLPF83vy4fvXZsfIaCYmatNxAACAOUaR5Q0ZHGjl0Ilz+eYPjjUdBQAAmGMUWd6QjzywNvN7urJ9t/FiAADgxlJkeUOWLZyXD9y7OjuGD6RW48UAAMCNo8jyhg0O9OXF42fz+P7jTUcBAADmEEWWN+zh9b3p6SrZPnyg6SgAAMAcosjyhq1YPC/vvWd1tu8eNV4MAADcMIosb8rWgVZ+cPR09hx4uekoAADAHKHI8qY8vKE3XSVWLwYAAG4YRZY3ZdXSBfnRO1f5niwAAHDDKLK8aVs3tvL0oVN5auxE01EAAIA5QJHlTdvS30opyTbjxQAAwA2gyPKmrV2+MO+6/WbjxQAAwA2hyDIthjb25YnRE/n+4VNNRwEAAGY5RZZpMTjQShJPZQEAgOtOkWVa3HLTorzttpuyY9j3ZAEAgOtLkWXaDA208t39x/P80dNNRwEAAGYxRZZpM9QeL9454qksAABw/SiyTJs7Vi3Jhr7l2bbb92QBAIDrR5FlWg0NtPKtH7yU0eNnm44CAADMUoos02poY18S48UAAMD1o8gyre5ZuzT3rl1qvBgAALhuFFmm3dBAK9949mgOnzzXdBQAAGAWUmSZdkMb+zJRk0dGxpqOAgAAzEKKLNPugdayrFu1ONuHjRcDAADTT5Fl2pVSMrSxL19++kiOnTrfdBwAAGCWUWS5LoYGWhmfqHl0r/FiAABgeimyXBcbb1mRW25alB3DtuEBAACmlyLLdVFKydBAK1986nBePnuh6TgAAMAsoshy3QxtbOX8+EQ+s/dg01EAAIBZRJHlunnHbTend/kCqxcDAADTSpHluunqKhnsb+WxJw/l1LmLTccBAABmCUWW62pwoC/nLk7ksScPNR0FAACYJRRZrqsH71yZVUvmGy8GAACmjSLLddXdVbK5v5XPPHEwZy+MNx0HAACYBRRZrruhgVZOnx/P579nvBgAAHjzFFmuu/fcvSorFs3L9uHRpqMAAACzgCLLdTevuysPb+jNp/aO5fzFiabjAAAAHU6R5YbYurGVE2cv5ktPH246CgAA0OEUWW6IH7tndZYt6Mn23VYvBgAA3hxFlhtiQU93PrJ+bR7dM5aL48aLAQCAN06R5YYZGujLsdMX8rXvH206CgAA0MEUWW6YD963JovmdWeb8WIAAOBNmFKRLaUMllKeLKXsK6X86uuc95OllFpK2dR+v66UcqaU8p32659ecu5j7Wu+8tnaN387zGSL5nfnww+syc6RsYxP1KbjAAAAHeqqRbaU0p3kd5MMJdmQ5KdLKRuucN6yJL+Y5GuXffR0rfXt7dfPX/bZxy/57OAbuwU6ydBAXw6fPJddzxovBgAA3pipPJF9MMm+WusztdbzSf4wyY9f4bzfTPJbSc5OYz5mmQ8/sDbze7qyfXi06SgAAECHmkqRvSXJ85e8398+9qpSyjuS3FZr/ZMr/PydpZRvl1I+V0p5/2Wf/cv2WPHfLaWUK/3yUsrPlVJ2lVJ2HTp0aApxmcmWLujJB+9bk50jo5kwXgwAALwBUymyVyqYrzaQUkpXkt9J8stXOO9Akttrre9I8ktJ/qCUsrz92cdrrRuTvL/9+pkr/fJa6+/VWjfVWjetWbNmCnGZ6YYGWjlw/Gy+s/+lpqMAAAAdaCpFdn+S2y55f2uSFy95vyzJQJLHSinPJnl3kk+WUjbVWs/VWo8kSa31m0meTnJf+/0L7T9PJPmDTI4wMwc8tL4387pLdhgvBgAA3oCpFNlvJLm3lHJnKWV+kp9K8slXPqy1Hq+1rq7MNniNAAAgAElEQVS1rqu1rkvy1SQfq7XuKqWsaS8WlVLKXUnuTfJMKaWnlLK6fXxekj+bZHha74wZa8Wiefmxe1Zn2+4DqdV4MQAAcG2uWmRrrReT/EKSnUn2JvmjWutIKeU3Sikfu8qPfyDJd0spjyf54yQ/X2s9mmRBkp2llO8m+U6SF5L8szdxH3SYoYFW9h87k5EXX246CgAA0GFKJz0R27RpU921a1fTMZgGR0+dz4/8T5/Kz3/wrvyNLQ80HQcAAJgBSinfrLVuutp5Uxkthmm3csn8vPuuldm+e9R4MQAAcE0UWRozONCXZw6fyvfGTjYdBQAA6CCKLI3Z0t+bUpJtuw80HQUAAOggiiyNWbtsYX7kjpW24QEAAK6JIkujhja28uTYiTx9yHgxAAAwNYosjRocaCWJp7IAAMCUKbI0qm/Forz9tpuyfdj3ZAEAgKlRZGnc1o2tDL/wcp4/errpKAAAQAdQZGnc0EBfkngqCwAATIkiS+NuW7k4/W9Znu2+JwsAAEyBIsuMsHVjX779g5dy4PiZpqMAAAAznCLLjGD1YgAAYKoUWWaEu9cszf29y4wXAwAAV6XIMmMMDrTyjWeP5uCJs01HAQAAZjBFlhljaGMrtSaPjIw1HQUAAJjBFFlmjPt7l+Wu1UtswwMAALwuRZYZo5SSwYFWvvrM0Rw7db7pOAAAwAylyDKjbN3Yl/GJmkf3GC8GAACuTJFlRul/y/LcevOibDNeDAAAvAZFlhmllJKhgVa+tO9wjp+50HQcAABgBlJkmXGGNvblwnjNp/caLwYAAP40RZYZ5+233pTW8oXZPjzadBQAAGAGUmSZcbq6Jlcv/tz3DuXkuYtNxwEAAGYYRZYZaWiglfMXJ/LZJw42HQUAAJhhFFlmpE3rVmb10vnZYbwYAAC4jCLLjNTdVbKlv5XPPnkwZ86PNx0HAACYQRRZZqyhgb6cPj+ez33vUNNRAACAGUSRZcb60btW5qbF87Jj+EDTUQAAgBlEkWXGmtfdlc0bevPpvQdz7qLxYgAAYJIiy4w2NNCXE+cu5kv7DjcdBQAAmCEUWWa0996zKssW9mT7bqsXAwAAkxRZZrQFPd356PrePLJnLBfGJ5qOAwAAzACKLDPe4EArx89cyFefOdJ0FAAAYAZQZJnxPnjfmiye351txosBAIAosnSAhfO68+EH1ubRPaMZn6hNxwEAABqmyNIRtg705fDJ8/nGs0ebjgIAADRMkaUjfOj+NVnQ05Xtuw80HQUAAGiYIktHWLKgJx+8b012jIxmwngxAADMaYosHWPrxr6MvXwu337+WNNRAACABimydIyPrF+bed0l261eDAAAc5oiS8dYvnBe3nfP6mwfHk2txosBAGCuUmTpKEMb+/LCS2ey+4XjTUcBAAAaosjSUR5e35vurpLtw8aLAQBgrlJk6Sg3L5mf9969Ktt3HzBeDAAAc5QiS8cZHGjl2SOn88ToiaajAAAADVBk6TibN7RSSowXAwDAHKXI0nHWLFuQB9etzPbdB5qOAgAANECRpSMNDbTy1MGT2XfwZNNRAACAG0yRpSMNDvQlSXYMeyoLAABzjSJLR2qtWJh33n5Ttu32PVkAAJhrFFk61tBAX/YceDnPHTnVdBQAAOAGUmTpWIMDrSRWLwYAgLlGkaVj3bZycTbeskKRBQCAOUaRpaMNbWzl8edfygsvnWk6CgAAcIMosnS0oVdXL/ZUFgAA5gpFlo525+oleaC1zDY8AAAwhyiydLyhgb7seu5YDr58tukoAADADaDI0vGGNrZSa7JzxHgxAADMBYosHe/etUtz15olVi8GAIA5QpGl45VSsnWgL1995kiOnDzXdBwAAOA6U2SZFQYHWpmoyaN7xpqOAgAAXGeKLLNC/1uW5/aVi7PNeDEAAMx6iiyzQiklQwOtfHnf4Rw/faHpOAAAwHWkyDJrDA60cnGi5lN7jRcDAMBspsgya7z9tpvylhULs334QNNRAACA60iRZdYopWTLQCuff+pwTp672HQcAADgOlFkmVW2buzL+YsT+cwTB5uOAgAAXCeKLLPKu26/OWuWLcj23caLAQBgtlJkmVW6ukq29PfmsScP5fR548UAADAbKbLMOlsH+nLmwng+9+ShpqMAAADXgSLLrPPgnStz8+J52T482nQUAADgOlBkmXV6uruyeUMrn3niYM5eGG86DgAAMM0UWWaloY2tnDx3MV986nDTUQAAgGmmyDIrvffu1Vm2sMd4MQAAzEKKLLPS/J6uPLyhN4/uGc35ixNNxwEAAKaRIsusNTTQl5fPXsxXnjnSdBQAAGAaKbLMWu+/d3WWzO/OjuEDTUcBAACmkSLLrLVwXnc+sr43O0fGcnHceDEAAMwWiiyz2tBAK0dPnc/Xnz3adBQAAGCaKLLMah+6f00WzuvK9t1WLwYAgNlCkWVWWzy/Jx+6b212joxmYqI2HQcAAJgGiiyz3tDGVg6eOJdv/eBY01EAAIBpoMgy633kgbWZ392VbcaLAQBgVphSkS2lDJZSniyl7Cul/OrrnPeTpZRaStnUfr+ulHKmlPKd9uufXnLuu0opu9vX/EellPLmbwf+tGUL5+X9967OzpHR1Gq8GAAAOt1Vi2wppTvJ7yYZSrIhyU+XUjZc4bxlSX4xydcu++jpWuvb26+fv+T4P0nyc0nubb8G39gtwNUNbezLCy+dyXf3H286CgAA8CZN5Ynsg0n21VqfqbWeT/KHSX78Cuf9ZpLfSnL2ahcspfQlWV5r/UqdfET2+0l+Yuqx4do8vL43PV0l24YPNB0FAAB4k6ZSZG9J8vwl7/e3j72qlPKOJLfVWv/kCj9/Zynl26WUz5VS3n/JNfe/3jVhOq1YPC/vuXtVdgwbLwYAgE43lSJ7pe+uvtoESildSX4nyS9f4bwDSW6vtb4jyS8l+YNSyvKrXfOHfnkpP1dK2VVK2XXo0KEpxIUr27qxL88dOZ09B15uOgoAAPAmTKXI7k9y2yXvb03y4iXvlyUZSPJYKeXZJO9O8slSyqZa67la65EkqbV+M8nTSe5rX/PW17nmq2qtv1dr3VRr3bRmzZqp3RVcweYNvekqyY5hqxcDAEAnm0qR/UaSe0spd5ZS5if5qSSffOXDWuvxWuvqWuu6Wuu6JF9N8rFa665Sypr2YlEppdyVyUWdnqm1HkhyopTy7vZqxT+b5D9M763BD1u1dEEevHNltiuyAADQ0a5aZGutF5P8QpKdSfYm+aNa60gp5TdKKR+7yo9/IMl3SymPJ/njJD9faz3a/uyvJvnnSfZl8knt9jd4DzBlWzf2Zd/Bk3lq7ETTUQAAgDeodNLCN5s2baq7du1qOgYdbOzls/nRv/fp/NLD9+UXH7q36TgAAMAlSinfrLVuutp5Uxkthlmjd/nCbLrj5mzbbRseAADoVIosc87gQCtPjJ7Is4dPNR0FAAB4AxRZ5pzBgVaSWPQJAAA6lCLLnHPrzYvztltXZPuw8WIAAOhEiixz0uBAX767/3j2HzvddBQAAOAaKbLMSUPt8eIdxosBAKDjKLLMSetWL8n6vuW+JwsAAB1IkWXOGhpo5ZvPHcvo8bNNRwEAAK6BIsuctXXj5HjxzhFPZQEAoJMossxZ96xdlnvWLrV6MQAAdBhFljlt60ArX//+0Rw+ea7pKAAAwBQpssxpgwN9majJIyNjTUcBAACmSJFlTlvftyx3rFpsvBgAADqIIsucVkrJ0EBfvvL0kbx0+nzTcQAAgClQZJnzhgZauThR8+ge48UAANAJFFnmvLfeuiK33LQoO4ZtwwMAAJ1AkWXOK6VkcKCVLzx1OCfOXmg6DgAAcBWKLGRyvPj8+EQ+88TBpqMAAABXochCknfefnPWLluQbbutXgwAADOdIgtJuromx4s/971DOX3+YtNxAACA16HIQtvgQCtnL0zksScPNR0FAAB4HYostD24bmVWLZlvvBgAAGY4RRbaerq7srm/N5994mDOXhhvOg4AAPAaFFm4xNBAX06dH8/nv2e8GAAAZipFFi7xnrtXZcWiedkxPNp0FAAA4DUosnCJed1d+ej63jy6dyznL040HQcAALgCRRYus3VjKyfOXsyXnj7cdBQAAOAKFFm4zPvuXZ2lC3qyY7fxYgAAmIl6mg4AM82Cnu48tH5t/v13XsiF8Yls7m/lg/etyaL53U1HAwAAosjCFf3K4APp7ir59N6D+XfffiEL53XlA/euyeBAKw890JsVi+c1HREAAOYsRRau4C03Lcpv/4W358L4RL7+/aPZMTyaR/aM5pE9Y+npKnn3Xauypb83m/tb6V2+sOm4AAAwp5Raa9MZpmzTpk11165dTcdgjpqYqHl8/0vZOTKWR0ZG88zhU0mSt992U7b0t7Klvzd3rVnacEoAAOhcpZRv1lo3XfU8RRauXa01+w6ezI7h0ezcM5rhF15Okty7dmkGB1rZ0t9K/1uWp5TScFIAAOgciizcQPuPnc4jI2PZOTKabzx7NBM1ueWmRdnc35st/a38yLqV6e5SagEA4PUostCQIyfP5dN7D2bnyGi+sO9wzl+cyMol8/PR9Wuzpb+VH7tndRbOswIyAABcTpGFGeDkuYv53JOHsmNkNJ994mBOnruYJfO786EHJkvth+9fk2ULrYAMAADJ1IusVYvhOlq6oCd/5q19+TNv7cu5i+P58tNH8sjIaB7dM5b/+N0Dmd/dlffesypb+lt5eENvVi9d0HRkAACY8TyRhQaMT9R86wfHsrO9WNTzR8+klGTTHTe3V0Bu5baVi5uOCQAAN5TRYugQtdbsPXAiO0ZG88jIaJ4YPZEk2dC3fLLUDvTm/t5lVkAGAGDWU2ShQz135FR2joxm58hYvvWDY6k1uWPV4lf3qn3HbTenywrIAADMQooszAIHT5zNo3vGsnNkLF95+nAujNesWbYgmzdMbuvz7rtWZX5PV9MxAQBgWiiyMMscP3Mhn31iclufx548lDMXxrNsYU8eaq+A/MH712TxfOu3AQDQuRRZmMXOXhjPF546nJ0jo/nU3rG8dPpCFvR05f33rsmW/t58dH1vbl4yv+mYAABwTWy/A7PYwnndeXhDbx7e0JuL4xP5+rNH88jI2KvFtrur5EfvXJkt/a1s7u9N34pFTUcGAIBp44kszCK11nx3//H2YlGjefrQqSTJ225dkc3tbX3uWbu04ZQAAHBlRouB7Dt4Mjvb2/o8vv94kuSetUuzpX9ysaiNt6ywrQ8AADOGIgv8kBdfOpNH2tv6fP3ZoxmfqHnLioXZ3B4/fnDdyvR0WwEZAIDmKLLAazp26nw+tXdyW58vPHUo5y5O5ObF8/LQ+sknte+/d3UWzutuOiYAAHOMIgtMyalzF/P57x3KzpHRfPqJgzlx9mIWz+/Oh+5fky39rXz4gbVZvnBe0zEBAJgDrFoMTMmSBT0Z2tiXoY19OX9xIl955kh2jozm0T1j2bZ7NPO6S95z9+ps6Z9cJXntsoVNRwYAYI7zRBa4oomJmm8/fyw729v6PHfkdEpJ3nn7za8uFnXHqiVNxwQAYBYxWgxMm1prnhw7kZ3Dk6V2z4GXkyQPtJZlS3tbn/V9y6yADADAm6LIAtfN80dPv7pX7a7njqXW5LaVi7JlQytbBlp55+03p7tLqQUA4NoossANcejEufYKyKP50r7DuTBes3rp/Dy8oZUt/b15792rM7/Htj4AAFydIgvccCfOXshnn5xcAfmxJw7m1PnxLFvQkw8/sDZb+lv50P1rsmSBNeYAALgyRRZo1NkL4/nSvsPZOTKaT+09mKOnzmd+T1fef8/qbOlv5aMberNyyfymYwIAMIMossCMcXF8IrueO5adI6N5ZGQsL7x0Jl0lefDOldnS38rm/lZuuWlR0zEBAGiYIgvMSLXWjLz4cnaOjGbH8GieOngySbLxlhWvbutzz9qlVkAGAJiDFFmgIzxz6OSre9V+5/mXkiR3rV6Szf2tDA608tZbVqTLCsgAAHOCIgt0nNHjZ/PontHsHBnLV545kvGJmtbyhdncflL74J0rM6/bCsgAALOVIgt0tJdOn8+n9x7MzpHRfP6pQzl7YSIrFs3LQ+snV0D+wL1rsmh+d9MxAQCYRoosMGucOT+ez33vUB4ZGc2n9o7l5bMXs2hedz5w3+QKyA890JsVi+c1HRMAgDdpqkXWho7AjLdofncGBya/M3thfCJfe+Zodo6Mtl9j6ekqec/dq7K5v5UtG3qzdvnCpiMDAHAdeSILdKyJiZrv7H/p1W19vn/4VJLkHbfflC39rWzpb+XO1UsaTgkAwFQZLQbmlFprnjp4MjuHR7Nzz2iGX3g5SXJ/77Js6e/N5v5W+t+y3LY+AAAzmCILzGn7j53OIyNj2TEyml3PHs1ETW69eVG2tLf1eeftN6fbtj4AADOKIgvQduTkuXxq71h2jozli08dzvnxiaxeOj8Pb2hlS39v3nv36szvsa0PAEDTFFmAKzhx9kIee/JQdoyM5rEnDubU+fEsW9iTjzywNoP9rXzw/jVZPN86eAAATVBkAa7i7IXxfGnf4ewcGc2je8Zy7PSFLOjpygfuW5Mt/a18dP3a3LR4ftMxAQDmDNvvAFzFwnndeWh9bx5a35uL4xP5xrPHXt3W59E9Y+nuKnn3XSsz2N/K5v5Wem3rAwAwI3giC3CZWmu+u/94doyMZufwaJ65bFufwf5W1tnWBwBg2hktBpgGtdbsO3gyO0dGs2PkP23r80BrWTa3S+36vmW29QEAmAaKLMB18PzR03lkz1h2jozmG88eTa3JbSsXZbC/lS39k9v6dNnWBwDgDVFkAa6zwyfP5VN7Jveq/dK+w7kwXrNm2YI8vKE3g/2tvPuuVbb1AQC4BooswA308tkL+ewTB/PIyFg+++TBnG5v6/PR9b3Z0t+bD9xnWx8AgKtRZAEacvbCeL741OHsGBnNp/aO5aXTF7JwXlc+cO+aDA608tADvVmxeF7TMQEAZhzb7wA0ZOG87nx0Q28+umFyW5+vf/9odoyM5pGRsTyyZyw9XSXvuXtVNve3smVDb9ba1gcA4Jp4Igtwg0xM1Dy+/6XsHJlcLOr7h0+llOQdt92UwYHJxaLuWGVbHwBg7jJaDDCD1Vrz1MGT2TE8mp0joxl58T9t6/NKqX2gZVsfAGBuUWQBOsjzR09n58hkqd313LHUmtyxanG2tLf1ecdtN9nWBwCY9RRZgA516MS5PNreq/bLT09u67N22YJs7u/N/9/evcfIdZZ3HP89c9+dvfm2u77mQkwu3lISQqCkSaEJtqEoQBVVUFGpVQVVBVVohdrCX1z+aREqVEWqmiZtUctFQIiKgNpOQ3NpqQOJE+r1JQ44ieN4d32Jd+29z+XpH+fM7szOrr3ZXfvMGX8/kjVnzr5n5lneGM9v3ve8745wW590km19AABA8yHIAkATGJko6LHnT2pX/6Aee/6UJgoldbakddcN3drR16s7t65TSyYZdZkAAAArgiALAE1mslDSE0dOadeBQT166KRGJgpqSSf1G28MtvV51w3d6mxhWx8AABBfbL8DAE0ml05q+7Zebd/Wq0KprKeOvjZzX+2uA4NKJUzvuG6tdmzr0btv6lF3O9v6AACA5rSoEVkz2ynpbyUlJT3g7n+1QLt7JX1H0lvd/emq81skHZT0WXf/UnjuJUnnJZUkFReTuhmRBYB65bLruePD2t0fBNqXz4zLTHrLllUzKyBvXt0adZkAAAAXtWJTi80sKemIpHdLOi7pZ5I+7O4H57Rrl/RDSRlJn5gTZB+SVJb01Jwge6u7n17sL0WQBYALc3c9P3Reu/uHtOvAoA4NBNv63LS+Qzu29WpnX6/e2NPGtj4AAKAhreTU4tsk/cLdj4Yv/C1J71cwwlrtC5K+KOlTcwr5gKSjksYW8V4AgGUwM93Q26Ebejt0391bdezM+MzU4688ekRf/s8junpNq3aEI7Vv3sS2PgAAIH4WE2Q3Snql6vlxSW+rbmBmN0va7O4/MLNPVZ3PS/oLBaO5NQFXkkvaY2Yu6R/c/f753tzMPibpY5K0ZcuWRZQLAKjYsqZVH73zWn30zmt18tyk9oTb+jz45Iv6h8ePqqcjq+03BSO1t12zmm19AABALCwmyM73Vf3MfGQzS0j6sqTfn6fd5yR92d1H55nGdru7nzCzbkmPmNlhd3+i7o2CgHu/FEwtXkS9AIB5dHfk9JG3X6WPvP0qjYwX9OPnh7Srf1DfeeYV/evel9XZktbdN/ZoZ1+v7ti6Vrk02/oAAIDGtJgge1zS5qrnmySdqHreLqlP0mNhWO2V9H0zu0fByO29ZvZFSV2SymY26e5fdfcTkuTuJ83sYQVTmOuCLABg5XW2pvXBmzfpgzdv0sR0SY8fOaU9Bwb1yMFBPbTvuFozSb3z+nXasS3Y1qcjx7Y+AACgcSwmyP5M0lYzu0bSq5I+JOl3Kz909xFJayvPzewxSZ8KF3u6o+r8ZyWNuvtXwynHCXc/Hx5vl/T55f86AIDXqyWT1M6+YHpxoVTW3qNntKt/UHsODulH+weVTpre8Ya12tnXq7tv7NG69mzUJQMAgCvcRYOsuxfN7BOSdivYfuef3P2AmX1e0tPu/v0lvG+PpIfDEdyUpG+4+64lvA4AYAWlkwndsXWd7ti6Tl94f5+efeWsdvUPaveBIX36e/v1Gduvt161Wtu39bCtDwAAiMyi9pFtFGy/AwDRcHcdGjiv3QcGtfvAoA4PnpckbdvQoZ3berWjr1dbu9nWBwAALM+K7SPbSAiyANAYXjo9NhNq9x0bliRduzav7eFetW/a2Mm2PgAA4HUjyAIALouhyrY+/YPae/SMimVXb0dOO8Lpx7dds1optvUBAACLQJAFAFx2w+PTevTQSe0+MKjHj5zSVLGsVa1p3XVjj3Zu69Wvs60PAAC4AIIsACBS49NFPXHklHb1D+rRwyd1frKo1kxS77q+Wzv6evWu69epnW19AABAlcUG2cVsvwMAwOvWmklpZ9967exbr+liWf8bbuvzyMEh/XD/gDLJhG6/bo12bOvV3Tf1aG0b2/oAAIDFYUQWAHBZlcqufcfOanf/oHYdGNTxsxNKmHTr1atnVkDe2NUSdZkAACACTC0GADQ8d9fBgXPaHe5V+/xQsK3Pr2zs1I5tPdrZ16vrutsjrhIAAFwuBFkAQOy8GG7rs6t/UM+9Em7rsy4fjNRu69WbNnWyVy0AAE2MIAsAiLXBkUntORjsVbv36GsqlV0bOnPavq1Xt1+3Vhu7WrShK6fOljThFgCAJkGQBQA0jbNj03r08Ent6h/Uky8E2/pUtKSTWt+V04bOFq3vzGl9V4s2zHlsy7K2IQAAccCqxQCAprEqn9G9b9mke9+ySWNTRT0/dF4Dw5MaGJnQicrjyKSOHDmlU6NTmvsdbXsuFQTdrlwQdsPQu6Fr9pH9bQEAiA+CLAAgVvLZlG7ZskraMv/Pp4tlDZ2b1MDInKA7PKnBcxPaf3xEZ8am665b1ZrW+s5guvL6MPRuqAq8PR05ZVKJS/zbAQCAxSDIAgCaSiaV0ObVrdq8unXBNpOFkgZHJnViZGJ2ZHdkUgPDEzp+dkI/ffE1nZss1lxjJq1tywbTleeM7lbCb3d7VqkkYRcAgEuNIAsAuOLk0kldvTavq9fmF2wzNlWsG9EdGJnQwMikfnFqVE++cEpj06Waa5IJU3d7tvZe3TmjvGvzWSUSLE4FAMByEGQBAJhHPpvSdd3tC+5j6+46NxmE3YHh2dHdyuOBV0f0yMEhTVctTCVJ6aSptxJwqwJvb9U05lWtrMQMAMCFEGQBAFgCM1NnS1qdLWnd0Nsxbxt312tj0xoYmdSJ4WA0t3o6889eOquhcwMqlmtXp8qlEzMLUs17325XTh259OX4NQEAaEgEWQAALhEz05q2rNa0ZdW3sXPeNuWy6/To1Mw9upXHSuj9n1+c1snzk5qTddWWTdVNYa4E3d7OnDZ05dSa4Z95AEBz4l84AAAilEiYujty6u7I6c2bu+ZtUyiVdfL8VH3QDR8Pnjin06NTddd1tqTrthmqHuXt7cwpm2LbIQBA/BBkAQBocOlkQhu7WrSxq2XBNlPFkoZGpoKpy9WLUw1P6sTIpPYdO6vh8ULddWvbMnX76q6fCbw59XTklGYlZiBWymXXeKGksali+Kek0amixqeLGg2fzx4XNTZdmmkXPJ+9ruSudMKUTiWUSpjSyYRSyeAxnag6TppSiYTSqYTSCVMqaUolE8okg+uC4+AxlTSlE+E14bXB6wbXVr9HKnzvuveqXFdzPqEki+ldMQiyAAA0gWwqqS1rWrVlzcLbDo1PF4P9dYfrtx568fSYfvLLMxqdqt12KGHSuvZs7b26c0Z517Zl+fAILEOxVA5C5HRtsJwNn8Hz8angeP4QOns8PmdF9QtpSSeVzyaVz6aUz6SUzya1Op/R5lWtymeTSiZMhZKrWCqrUHYVimUVy65CqaxCqazpYllj06XwfFnFkqtQLqtQdBXLZRVKQdvKefeL17QcZqoLyUHArg29Fw7NlWurrk/OCfF15y/eNpW0INiH75VJVYf84OephLHY3yIRZAEAuEK0ZlJ6w7o2vWFd24Jtzk8WaqYtz0xnHpnQ4YHz+vHhk5os1K7EnEqYejpy9YtShUG3tzOnNfkMH87QFNxd05XgOWf0cjZ8lmaCZuW4OpDOXhc8n5qzuvlCzDQTNoPH4Li3IzdzXH0+n02pLZtSa3hN5bgtm1Jr2PZyfgnl7iqVfSYIF8OQWyiHQblUDkNz8L9xsVQdmueG6apzYUieeb1S1euFrz37mkHAni7Nnp8olGZD+Nz3qqqzWA7qv9TqRr6rAnJtwJ4diV44NAfPK6H5hvUduudXN1zy3+FyIMgCABS9JYgAAAkJSURBVIAZ7bm02nNpvbFn4W2HhscLdSO6lcD77Ctn9R/9kyqUaj/sZVKJmenKlZWXK6O8q1ozMx+8KqMTqXBqYjIRjKIkwymEyYQpnWTEAovn7poslKtGMOun0I5OlYLRzumixi8SSMeminUrjS8kmTDlM0GAzGdTas2m1JZNanW+NTxXHTxTymeqw2dy9nwYQnOpZKz3oTarTDkO9vOOo3J5ntAcjkDXnp8nhM8E4tnAXtO2OE+oL88TwoteE7zHp4vhe1WCd/V71YbxHX29BFkAAHDlMTOtyme0Kp/Rtg0Lr8R8Zmy69l7dqlHevUfPaOj81LJGNhKm2cBbHX7rgnAwQjETiCv37tW0q742aJtKJOpea27bZDg1MRmOilRCdrI6lC90TFhfULnsGpuuDY61wbN+2u34TCCdEz7D6bqL/U8tk0rMhMjKY3suWCE8GMlM1oXOC4XQbCpxxfZjs0okTNlEUtmYpii/1HO7L6OYdgEAAGhUiYRpXXtW69qzetOm+dsUS2WdGp3SieFJnZsoqFiZ/leeHWkohlMQi6WySuFoQykcyShV/axyzWybYOSh+pqgbdBusuizr19pVy6rVHIVyrXXV665DLMJ68wX1oPQazWhty5sVwXmmiBfE7ZXPqxXAnmxPOd+z3kWEwpGPucuQhT8bKLweu/vDAJmZcrs6nxGm1e3qi0TTJ9tmxM0K+3yVaG00paFzdDsmumLFYIsAAC47FLJRLhw1MIrMTeSciU4h6G4FN6TNxt2q0J1abbd3OPqID4Ttsuu0kwgn/s69WF9bsife02p7JoslmqDfN1r1R9f6rBuppnAWL2wUGXP47mBtPqez9pAGpxvvcz3dwJoLARZAACAi0gkTJmEKaPmHbGrhPXaEera8DxvEC/VXpNMWDjymawJpC3pZFONBgGIFkEWAAAAM2FdkloUz4V4AFw5mvdrRQAAAABAUyLIAgAAAABihSALAAAAAIgVgiwAAAAAIFYIsgAAAACAWCHIAgAAAABihSALAAAAAIgVgiwAAAAAIFYIsgAAAACAWCHIAgAAAABihSALAAAAAIgVgiwAAAAAIFYIsgAAAACAWCHIAgAAAABihSALAAAAAIgVgiwAAAAAIFYIsgAAAACAWCHIAgAAAABihSALAAAAAIgVgiwAAAAAIFbM3aOuYdHM7JSkl6Ou4wLWSjoddRFYMfRnc6E/mwv92Vzoz+ZCfzYX+rO5xKE/r3L3dRdrFKsg2+jM7Gl3vzXqOrAy6M/mQn82F/qzudCfzYX+bC70Z3Nppv5kajEAAAAAIFYIsgAAAACAWCHIrqz7oy4AK4r+bC70Z3OhP5sL/dlc6M/mQn82l6bpT+6RBQAAAADECiOyAAAAAIBYIciuADPbbGb/ZWaHzOyAmd0XdU1YOjPLmdlPzeznYX9+LuqasDxmljSzZ83sB1HXguUzs5fMbL+ZPWdmT0ddD5bHzLrM7Ltmdjj8d/TXoq4JS2Nm14d/Lyt/zpnZJ6OuC0tjZn8afg7qN7Nvmlku6pqwdGZ2X9iXB5rl7yVTi1eAma2XtN7d95lZu6RnJH3A3Q9GXBqWwMxMUt7dR80sLem/Jd3n7nsjLg1LZGZ/JulWSR3u/r6o68HymNlLkm5190bfBw+LYGZfk/Skuz9gZhlJre4+HHVdWB4zS0p6VdLb3P3lqOvB62NmGxV8/rnJ3SfM7NuSfuTu/xJtZVgKM+uT9C1Jt0malrRL0h+7+wuRFrZMjMiuAHcfcPd94fF5SYckbYy2KiyVB0bDp+nwD9/4xJSZbZL0W5IeiLoWALXMrEPSnZIelCR3nybENo27JP2SEBtrKUktZpaS1CrpRMT1YOlulLTX3cfdvSjpcUkfjLimZSPIrjAzu1rSzZKeirYSLEc4FfU5SSclPeLu9Gd8fUXSn0sqR10IVoxL2mNmz5jZx6IuBstyraRTkv45nP7/gJnloy4KK+JDkr4ZdRFYGnd/VdKXJB2TNCBpxN33RFsVlqFf0p1mtsbMWiW9V9LmiGtaNoLsCjKzNkkPSfqku5+Luh4snbuX3P3NkjZJui2ckoGYMbP3STrp7s9EXQtW1O3ufouk90j6uJndGXVBWLKUpFsk/b273yxpTNJfRlsSliucIn6PpO9EXQuWxsxWSXq/pGskbZCUN7OPRFsVlsrdD0n6a0mPKJhW/HNJxUiLWgEE2RUS3kv5kKSvu/v3oq4HKyOc4vaYpJ0Rl4KluV3SPeE9ld+S9Jtm9m/RloTlcvcT4eNJSQ8ruOcH8XRc0vGqWS/fVRBsEW/vkbTP3YeiLgRLdrekF939lLsXJH1P0jsirgnL4O4Puvst7n6npNckxfr+WIkguyLCxYEelHTI3f8m6nqwPGa2zsy6wuMWBf9nfjjaqrAU7v5pd9/k7lcrmOb2Y3fnG+UYM7N8uKiewimo2xVMmUIMufugpFfM7Prw1F2SWCgx/j4sphXH3TFJbzez1vBz7l0K1oBBTJlZd/i4RdJvqwn+jqaiLqBJ3C7p9yTtD++rlKTPuPuPIqwJS7de0tfCFRcTkr7t7mzbAjSGHkkPB5+rlJL0DXffFW1JWKY/kfT1cDrqUUl/EHE9WIbw/rt3S/qjqGvB0rn7U2b2XUn7FExBfVbS/dFWhWV6yMzWSCpI+ri7n426oOVi+x0AAAAAQKwwtRgAAAAAECsEWQAAAABArBBkAQAAAACxQpAFAAAAAMQKQRYAAAAAECsEWQAAGoyZjVYdv9fMXgj3/gMAAGIfWQAAGpaZ3SXp7yRtd/djUdcDAECjIMgCANCAzOwOSf8o6b3u/suo6wEAoJGYu0ddAwAAqGJmBUnnJb3T3f8v6noAAGg03CMLAEDjKUj6iaQ/jLoQAAAaEUEWAIDGU5b0O5LeamafiboYAAAaDffIAgDQgNx93MzeJ+lJMxty9wejrgkAgEZBkAUAoEG5+2tmtlPSE2Z22t3/PeqaAABoBCz2BAAAAACIFe6RBQAAAADECkEWAAAAABArBFkAAAAAQKwQZAEAAAAAsUKQBQAAAADECkEWAAAAABArBFkAAAAAQKwQZAEAAAAAsfL/Qkhmpvjx9qQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x648 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "CVs = []\n",
    "for k in k_range:\n",
    "    f = open('admix.%d' % k)\n",
    "    for l in f:\n",
    "        if l.find('CV error') > -1:\n",
    "            CVs.append(float(l.rstrip().split(' ')[-1]))\n",
    "            break\n",
    "    f.close()\n",
    "fig = plt.figure(figsize=(16, 9))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(k_range, CVs)\n",
    "ax.set_title('Cross-Validation error')\n",
    "ax.set_xlabel('K')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load meta-data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = open('relationships_w_pops_121708.txt')\n",
    "pop_ind = defaultdict(list)\n",
    "f.readline()  # header\n",
    "for l in f:\n",
    "    toks = l.rstrip().split('\\t')\n",
    "    fam_id = toks[0]\n",
    "    ind_id = toks[1]\n",
    "    if (fam_id, ind_id) not in ind_order:\n",
    "        continue\n",
    "    mom = toks[2]\n",
    "    dad = toks[3]\n",
    "    if mom != '0' or dad != '0':\n",
    "        continue\n",
    "    pop = toks[-1]\n",
    "    pop_ind[pop].append((fam_id, ind_id))\n",
    "#ind_pop[('2469', 'NA20281')] = ind_pop[('2805', 'NA20281')]\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def load_Q(fname, ind_order):\n",
    "    ind_comps = {}\n",
    "    f = open(fname)\n",
    "    for i, l in enumerate(f):\n",
    "        comps = [float(x) for x in l.rstrip().split(' ')]\n",
    "        ind_comps[ind_order[i]] = comps\n",
    "    f.close()\n",
    "    return ind_comps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "comps = {}\n",
    "for k in k_range:\n",
    "    comps[k] = load_Q('hapmap10_auto_noofs_ld.%d.Q' % k, ind_order)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ordering = {}\n",
    "for k in k_range:\n",
    "    ordering[k] = cluster(comps[k], pop_ind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 648x648 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAHICAYAAABXk5v/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuQW2X9wOFvfgVBqAhYlKKwFVRAhdGptkVlmi0KFhXUUkRU8I4jI3fLKCJbGLmKOsiAwgxex1GsgGLHa+mig0OpHUZFvAIChXIRFC1YpW1+f5Rd9pLsJtk3OZc8zwzDNskmb86eJOeT9+SkUqvVAgAAgHT+L+sBAAAAlI3QAgAASExoAQAAJCa0AAAAEhNaAAAAiQktAACAxIQWAABAYkILAAAgMaEFAACQ2FatXLiyXaUWO275efa6TgynGNbE7IiZa8af8UREbPf0P4u8jNbMTHdds9dFPBIRz0l3laWzZmaMW39Ib/q/p8f6Z60ff8bIZb9u9tOn13ucj7xMo/PzrNHY1415Xhu5HJ6yzTa3xdZbbz3qtPXr947p0//U8jCa/b316/eOiGh42aHrGbrc0x6JmPm38fdj5pqIdbO3/M6z1j99v+vc36HbHXfdQ8upwe80beztDv273t+mYzw7t2bi5TW0no5fH+tftt7jYORpjdfvqHuZLf+fEdOn/33c5dt9rNLYyL/L038nj6kymz074pFHHom//e1vf6/VartMdvlKrVZr+soru1VqcdyWn2sD7Q6x+CpRixiojD9jZUTtxu6PpxMqA+muqzYQMRBb/qO+ykBErIyI/owHUnLVwWoMVgfHnzFy2Q+MeE6s9ziPEZeb6PyiGxj/2tDXNyte+MIXjjptcHBlVKutr7jN/t5klxs6f3BwZUTEiJ8HIgaWjr8fA5WJ/3ZjLj/yusf9Tp1l1JKhsUw2tqnezsSDCM/OrRiIestraD0ZWleH15kJjP2deia6zpGnj7yOwcH5Ua2WZGMk50b+TUY993hMlVatFjEwMBBLly5dU6vVXjXZ5e06SF29HNKUV93IaleZI2usgUry+9tOnOXOVJZJw9DrofWqB3VivR97nX19X01+G0B7hBYN1QZGB5f4gh7Uaxv+bdzfpjee6wXr0L+FbOF1I6KAYmnpM1qMNjI8Uu5qlzcCC3pQTgKr2Q3NljZIJ7tvY3fpa+Z3WtGlZTtqt0c6ot1lLKCgN5jRasdARXzE+NmusTNgUGg5CY08G6wOZr7BmPXtDxuarWpnvWkm/EaY8n22bmcqN+ss0HFCqw1i4mmTxZVlReF04PNIdFmj3fPa0Pbn+lKvQ9bL3OpmOIm0/PE3YSJCi44wuwVkav7SLf8fGSdZhkreIkm4dZWN8XLz96WRlkPLxjOTsY4AndbyTFM3o6KTtzXy4BkT3M6o5TPysuIKkhFYTKatGS0b0lgHgKlIeqj9Dki6AdWN8GrmfJGVqamuUzbqoXjsOkhL7BIIvevuWXd3/TbzHmST6kDc1FsmwxvhdgnMlcniSDxBubUUWrPXdWoYAFBMDWMw4+ApfKQCFJwZLQBKoyNxUZQZoqKMsyTMVgGTEVoAdJWZlrSq1X7LtEv6+r4aESIKaI7QSsTnlgA6aygmxv4f8kyUQe8SWgC0pBcCp+n7aHe9niOcgGYJrYTOujHrEQAAZSb0oDiEFgClkHymzWwVOSW2oBiEFgBT1ihyBquDo87rxG6HY2+j1/TyfQfIM6EFQEcIAAB6mdACoCtaCa9Oz4IBQKcJLQBaNlH8pAojgQVAkW2V9QAAKLaU32s10We9AKBIzGgBAAAkJrQAaFq7M0uT/V6nrhey4PDrQITQAiAnmokmYQVAUQgtAACAxIQWAABAYkILgLbZlQ8A6hNaAAAAiQktADJlVgyAMhJaACQnngDodUILAAAgMaEFAACQmNACAABITGgBkAmf4wKgzIQWAABAYkILgMyY1QKgrIQWAABAYkILAAAgMaEFAACQmNACAABITGgBAAAkJrQAAAASE1oAAACJCS0AAIDEhBYAAEBiQgsAACAxoQUAAJCY0AIAAEhMaAEAACQmtAAAABITWgAAAIkJLQAAgMSEFgAAQGJCCwAAIDGhBQAAkJjQAgAASExoAQAAJCa0AAAAEhNaAAAAiQktAACAxIQWAABAYkILAAAgMaEFAACQmNACAABITGgBAAAkJrQAAAASE1oAAACJCS0AAIDEhBYAAEBiQgsAACAxoQUAAJCY0AIAAEhMaAEAACQmtAAAABITWgAAAIkJLQAAgMSEFgAAQGJCCwB6VF/fV7MeAkBpCS0AAIDEhBYAAEBiQgsAACAxoQUAAJCY0AIAIiKiWu3PeggApSG0AAAAEhNaAMCwarXfzBZAAkILAAAgMaEFAD3M7BVAZwgtAGBccAkwgKkRWgDQ40QVQHpCCwCYkANkALROaAEATRFcAM0TWgAAAIkJLQAAgMSEFgDQkF0FAdojtACAutqJLGEGsIXQAgA6ykE0gF4ktACArhBbQC8RWgBAUhMFldktoFdslfUAAIDyEVNArzOjBQAAkJjQAgAASExoAQAAJOYzWgBAS0Z+/mpwcGWGIwHILzNaAEDbHEUQoD4zWgDAlIktgNGEFgBQGGODzq6LQF7ZdRAAKCwzaUBeCS0AoNB8TgzII6EFAACQmNACAHLD7BRQFkILAMglwQUUmaMOAgC5MxRZvhwZKCozWgBAIXRrhstMGpCCGS0AgDrGBtddd2U0EKCQzGgBAIxhVguYKqEFAPQ8YQWkZtdBAIAQW0BaZrQAgJ4msIBOEFoAAACJCS0AoGfU+34ugE7wGS0AoCeMjSyxBXSSGS0AAIDEhBYAUCpmqoA8EFoAQOmMjC3hBWRBaAEApSe2gG5zMAwAoJTEFZAlM1oAQGmIKyAvhBYAUAoiC8gToQUAAJCY0AIAus7sE1B2QgsAACAxoQUAZMbMFlBWQgsAyJTYAspIaAEAACQmtACATJjJAspMaAEAACQmtAAAABITWgBALtiVECgToQUAAJCY0AIAAEhMaAEAACQmtACAzPl8FlA2QgsAACAxoQUAAJCY0AIAAEhMaAEAACQmtAAAABITWgAAAIkJLQAAgMSEFgAAQGJCCwAAIDGhBQAAkJjQAgAASExoAQAAJCa0AAAAEhNaAAAAiQktAACAxIQWAABAYkILAAAgMaEFAACQmNACAABITGgBAAAkJrQAAAASE1oAAACJCS0AAIDEhBYAAEBiQgsAACAxoQUAAJCY0AIAAEhMaAEAACQmtAAAABITWgAAAIkJLQAAgMSEFgAAQGJCCwAAIDGhBQAAkJjQAgAASExoAQAAJCa0AAAAEhNaAAAAiQktAACAxIQWAABAYkILAAAgMaEFAACQmNACAABITGgBAAAkJrQAAAASE1oAAACJCS0AAIDEhBYAAEBiQgsAACAxoQUAAJCY0AIAAEhMaAEAACQmtAAAABITWgAAAIkJLQAAgMSEFgAAQGJCCwAAIDGhBQAAkJjQAgAAekq12t/x2xBaAAAAiW2V9QAAABjve9/7Xuy8885Tuo5HH300Fi1alGhEQCuEViPr1kXsumtLvzIw9oQHHoiYOTPViPLDsiHH1p26Lnad3tr6GWeN/ucD6x+ImRdbP/PC35ROKELETHV8qa4DyqqNTdoYGBiIgYGB2RFRe+qkByOi7rUIrUZaXeqduo48smzIsZY3yDt0HaTjb9qcFEHaS7MfIgZItDn6vEZnCC2AHmSWqHxSxKRwgKe1N+s5MOpfduDpbUILOqHFuehahI1YusosEcDEUrzxYAee3uaog9AJCZ5ZbcQCABSX0GrBhg0b4q677ooNGzZkPRQAIKH169fXPf1Pf/pTl0cCZGXTpk11T3/ooYfauj6h1YSHHnooDjvssHjjG98Yn/zkJ+OQQw6Jww47rO2FXjbr1q2Lf/zjHxERcdNNN8X1118fGzduzHhUQBmdeeaZ8Z///GfUaffee2+8853vzGhElMXChQvjkUceGXXaypUr46ijjspoRKO94hWviGOPPTYuv/zyuPXWW2Pz5s1ZDwlK561vfWv897//HXXa7bffHoccckhb1+czWk1YsmRJnHDCCfH6179++LQVK1bExz/+8fja176W4ciy9/GPfzxuvfXW2LhxY/T19cWznvWs2HHHHeOb3/xmfOc738l6eJm77LLL4qMf/WjWw+hZS5cujUqlMvzvZz7zmbHffvvFG9/4xgxHxVTsu+++sWDBgjjxxBPjsMMOi/PPPz9++ctfxrnnnpv10Ci4888/P970pjfFd7/73dh9993ja1/7WlxxxRXxox/9KOuhRUTEmjVr4re//W2sWrUqLrnkkrjtttti+vTpMXfu3Dj//POzHh6UwrHHHhuHHnpoXHvttbHDDjvEihUr4rTTTotvfetbbV2f0GrC2rVrR0VWRMRBBx3kiS0ifvWrX8VNN90UmzZtiv322y9uv/32iIjo7+/PeGT5sGzZMqGVoXnz5o3694YNG+JnP/tZ/PCHP4xLL700o1ExFUcffXQcfvjh8frXvz4+8IEPxKc+9alYuXJl1sPKpbFvNIz06U9/usujyb/Xvva1ceWVV8YRRxwR8+bNi3vvvTdWrFgR2267bdZDi4iIadOmxStf+cqYMWPG8H+/+c1v4g9/+EPWQ+sJl112WcPzvM6XxxFHHBE77bRTHHroobFo0aK49tpr46c//WnssssubV2f0GpCrVare7pp+4htttkmIra8AOw64gAQjV7ce83atWsbPjl7Yu68elP9hx9+uDcCxnj1q1896jE79JxXqVTilltuyWpYdd16662xZMmSWLBgQZx33nmxdOnSePLJJ2PJkiW52SDOi5FvNJx++ulx4YUXNnw9Y8seGpVKJXbffff4+te/Hh/4wAeGg/TCCy/MeHQRixYtikqlEnvttVfMmTMnTjjhhNh9992zHlbP2H777euebnunXIa22V70ohfFOeecE2eccUZ897vfjYj2ttuEVhPuvffemDNnzqiNj1qtFv/6178yHln21q5dO7xsHnvsseGfLZstpk2bFtOnT7dxkyMbNmyI//3vf1kPI1dWr16d9RCadu6558aVV14Zs2bNioiIarUaX/nKV2L+/PmxatWqbAeXMyPfaDj//PPj4IMPznA0+ffmN795+OePfexjGY6kvte85jXxm9/8ZvigXP/9739j7ty5sddee2U9tJ5w7LHHjjtt8+bN8b3vfS+D0dApQ0Hd398/6k3ZdoNaaDXhz3/+c9ZDyC3LZmK77rprHHPMMVkPo2ctXrx41JPjhg0bYu3atXH22WdnOKr8Wb58eRx44IGxww47xB133BFnnnlm1Gq1+PSnPx377rtv1sMbZeidxZHe9773xeLFizMYTXF4131yN910UxxwwAExZ86chrMXWTr11FOHf3744Yfj+uuvjyOOOCLuvffe+Pvf/57hyHrDrbfeGuedd15EbNn19uabb44rrrgi3vCGN3j+KZHdd989FixYMO70Cy64oK3rE1pN+PrXv97wvF7fiLZsJnb00UfH/fffH/fdd1/sueee8ZznPCfrIfWUz372s6P+ve2228bznve8jEaTX+ecc07cfPPNEbHlXdvPfOYzMWPGjDj++OPjhhtuyHh0ox122GHxgx/8ICIiPvzhD8cVV1wxfHrexpq1oTcaarVa3HbbbXHkkUdGrVaLSqUSV199ddbDy52+vr647rrr4owzzoiNGzfGK1/5ypg3b17MmzcvF284rF69OlatWhU333xz3H777bHTTjvFwoULY+7cuVkPrSccf/zxcfHFF8c//vGPWLBgQZx55pnxq1/9KrbayqZ0mZxzzjnx2GOPxdve9raI2DJredxxx8W0adPauj5rRxNGHkr485//fJxyyil2BXuKZTOx++67L9785jfHPvvsE7/73e/i5JNPjve///1ZD6tnLF++fHif6muuuSbe/va3R8SWz1ssWbIky6HlyjOf+cyI2PJVFo8//njMnz8/4xE19u9//3v457/85S/DP3veGW/sGw1M7F3vele8613viogts9/f+ta34pJLLokPfvCDDb9bp5u+/OUvx7x58+L000+Pl7/85WYpu2zbbbeNAw44ICIi9txzz1zuXsrULV++PN7xjnfEo48+GkceeWQsXrw4FixY0PY2g9BqwnHHHTf887e//e348Ic/nOFo8sWymdhPfvKTWLNmTVQqlfjPf/4TCxcuFFpdNPKoj5deeulwaP34xz8WWiNUKpW4/PLL45ZbbokjjzwyIiKefPLJcd9XlQcjNy4b/cwWjz/+eNZDKJTrr78+Vq1aFXfeeWdEROyzzz5x3nnnxZw5czIe2RZLly6NGTNmxDbbbBMbNmyIb3zjG7F58+Z4z3veE9ttt13Wwyu9kbPCd955pxniktpuu+3i2muvjWOPPTYuvPDCOPfcc2PRokVtX5/QapEX88Ysm/E2bdo06tC7TzzxxPAh8F/60pdmNSwY5bHHHotarRbz58+P97znPRGx5YvIP/WpT2U8svGGDsATEfHPf/5z+OfHHnssy2Hl0kUXXTS86+DY5+errroqo1Hl1xlnnBG77LJLLF68OObOnRv7779/27sLdcLixYuHv8rg/e9/f7zkJS+J5z73ufHud787rrnmmoxHV377779/fPKTn4yddtop66HQQUNH4d28eXNUKpW44IILhj+f1c5ReIVWE4YO+Vqr1eKvf/3rqHfC83DI1yxZNhN72cteNmr3nZe//OXD/7ah03lD6+TI9bNWq8Udd9yR9dBy5dnPfva4w9busccesccee2Q0osZOPPHErIdQGG95y1tit912i3nz5sXhhx8e999/f0REnHzyyRmPLJ9++9vfxuOPPx6rV6+On/3sZ/G5z30unnzyyXj+858fF198cdbDi2c84xmxzTbbxBNPPBGrV68e/gLVZcuWZTyy3nDFFVfEkiVLYsaMGXHuuefGzjvvnPWQaFO12vgrXt773vcmvS2h1YRZs2ZFrVaL/fbbb/jwr7fddluu3unKimUzsfXr19c9fdq0aXHKKafEhz70oVx8yLqsTj/99Lrr5yc+8YmMR5Yvv/71r8ftHjU0C5K379E688wz4yUveUm87W1vi5kzZ0ZE1J2xIeKLX/xi/PznP4+IiH/961+xevXq2LRpUxx66KFx9NFHZzy6fNq0aVNs3Lgx/ve//8U///nPuOeee3LzdRC1Wi3+8Ic/xA033BBveMMbhk9/4oknMhxV79hzzz1j2bJlceqpp8a+++4bfX19uX2epH2pX2OEVhOuu+66uPrqq0dNF++///5x1FFHxUc+8pEMR5Y9y2ZijT6MXqvV4u677473ve99w0d7Iz3rZ3Nmz549vEtS3j344IPx85//PK699tq4+eabo7+/P4444ojYbbfdsh5aLg296TX0xbvTpk2LjRs3Zjmk3Jo9e3Zsv/32MWfOnJg7d24cc8wxsccee+RiNitiy+dMzzrrrNh+++3joosuioiIO+64Y9T3f9E5f/zjH+O0006LvffeO/7yl7/EDjvskPWQ6IDUrzFCqwlPPvnkuH1yd9ppp9y8y5Uly2ZifX19Dc+bNWtWLF26tIuj6T3Wz/LZeuutY+HChbFw4cL4xS9+ESeffHLceeed8YUvfCHroeVOpVKJhx56KJ773OcOf/HmunXrYvPmzRmPLJ9WrVpV91DdP/nJT0Z9h1VW9ttvv7j66qvjzjvvjKuuuiqWL18ee++9d93v/CG90047LT772c/GPvvsk/VQ6KDUrzFCqwmVSiUefvjh2GWXXYZPe/DBBx1OOCybqTrkkEOyHkKpWT+b89Of/jTrITTtxhtvjGuuuSbuueeeOOCAA+Lqq6+OvfbaK+th5dI555wTBx98cLz97W+PmTNnxtq1a+O6666LL33pS1kPLZfy/n1IF1xwQaxYsSL23HPPOOqoo2LFihXD3yNH5/3whz/Megh0QerXmHw/q+TE2WefHQcddFAsWrQodtttt+EXq8suuyzroWXOsiHPrJ/N2XrrrbMeQtP6+/vjVa96VfT19cXq1avj17/+9fB5DrE82mtf+9q44YYbYvny5bFu3bp48YtfHCtXrvQh/gaGvuB5pFqtFr/73e8yGtFoP/jBD6Kvry8OP/zweN3rXuez0NABqV9jhFYTDjzwwBgcHBz1YjU4OOgQn2HZkG/Wz/K56667sh5Coey8887Dh+xnYnn/guebbrop1q5dG8uWLYvLLrssfv/738f3v//9mD9/fuy4445ZDw9KIfVrjNBqkherxiybOh54IGLXXad2FesfSDSY3mb9rO+B9Q/ErtOLt45O9LlHmIoirFsveMEL4qSTToqTTjop7rvvvli2bFl87nOfixtvvDHroZXSo48+OuUZ4Ae8lBdK6ucBoQWd8NQhQZtVGYiIlRHR+KsdIKmZF7e2jkZEnDUQsXQg+VCANjz/+c+PE0880XfLddCiRYtauvzg4MqIGHjqP4j4v6wHAJBSilkWs4kUUYr19tFHH00wEgAizGg1lmDXr9LOF1s25Fg7MzVmE/OtqLs5dls7637fV/rihS98YQdGk38pdgvrdJgWYYxQZCkeYxHxYKMzhFYjLe76FdFDk8WWDdBF4plOaHW3sIgtH5TvZpgWYYxQFNVq/1O7dz5t0aJF406bSK0WMTAwEEuXLl1Tq9VeNdnl7ToIQEtqA1mPAADyT2gBAAAkJrQAAAASE1oAAACJCS0AAKD0qtXuHiVJaAHQMgfEAICJCS0AAIDEhBYAAEBiQgsAACAxoQUAAJCY0AIA6LLB6mDWQwA6TGgB0DRHGwSA5ggtAACAxIQWAABAYkILAAAorWq1P5PbFVoAAACJCS0A2uLAGADQmNACAABITGgBAAAkJrQAmBK7EALAeEILAAAgMaEFAACQmNACAABITGgBAAAkJrQAAAASE1oAtM0RBwGgPqEFAJCBwepg1kMAOkhoAQAAJCa0AAAyZGYLykloAQAAJCa0AAAAEhNaAAAAiQktAJKpDTjkOwBECC0AAIDkhBYAAEBiW2U9AACKz+6CADCaGS0AAIDEhBYAQEZ8WTGUl9ACAABITGgBAAAkJrQAAAASE1oAAACJCS0AAIDEhBYAHeG7tQAoomq1P8n1CC0AAIDEhBYAQMZ8nxaUj9ACAABITGgBAAAkJrQAAABGSHFADKEFAACQmNACIDmHdgeg1wktAACAxIQWAABAYkILAAAgMaEFAACQmNACoGuKcpCMoowTgM6Z6iHehRYAAEAd1Wp/28EltADoOjNGAJSd0AIAAEhMaAGQuTzOcOVxTAAUh9ACIJeEDgBFJrQAAHJmsDqY9RCAKRJaAHRVMzNVZrMAKDqhBUCm8hxVeR4bAPkmtAAAABITWgAAOeRzWlBsQguAXLG73miWB0AxCS0AcqE20NtR0cv3HSDvqtX+ln9HaAFAAQgxgGIRWgB0TBnioJP3od51jzxt7PllWJ4AvUJoAZAZ4bBFvaCaKLgmO53ycEAMKC6hBQAZaiaoGl0egPwSWgDQZVONpV4/cAjQmnYO5NCLUi8noQVARwmC+rJaLv4exWP3QSimrbIeAAD0ik4fWKMy8PRtVCa4raHLzOrr3HgAiijlrJYZLQDogm4cvdBsVbmZ2YJiEVoAUFLCCyA7QgsASsyBM8rBbBZ0VicOGCK0AMhEETf+8zjmZr9jS3AB1NepozIKLQDoMIFDs8xcQXkILQA6rpdDo5fvO60RWVAuQgsAWjDREf5EFSkILigHoQUALRJZAExGaAHAFDV7QAoA8qVTB8KIiNiqY9cMAD1KYAFgRguArih6fLRyGHUAEFoA0CZRBUAjQgsAJiGoAGiV0AIAAEhMaAEAACQmtADoGXYBBKBbhBYAhdJKLAkrALIitAAoNV8mDEAWhBYAPUtsUVSD1cGshwBMQmgB0FPEFWUhtiDfhBYAhdNqLNUGRv+O2AJ6WbXan/UQeoLQAgAASExoAVBKZq0AyJLQAqCQhBSM5jNbkC9CCwCg4EQW5I/QAqBUzHQBTM4BMTpPaAEAQMkJq+4TWgAABWa3QZo1Nraq1f7h04RYekILgNI568asRwBArxNaAAAAiQktABjBwTQAyiPLXSKFFgCFliqMBBYAKQktAHiK2ALKzkEvukdoAQAUlCMOQn4JLQAKz0wUQPvKOMuVh/sktAAAgJaM/A4u6hNaAJSGmS0A8kJoAdDTxBlF4jNZdEOrM1VmtuoTWgAAAIkJLQAAgMSEFgAAEBF2A0xJaAFQClP5rJXPaQE0Jr7aI7QAKCyBBEBeCS0AAGDSmat655vtakxoAQC0waHW6UVFCqusxyq0AACmoNvBJfDIQtbRUkRCCwCgTaIHaERoAQBMkeCirMxktU9oAQAApZGXONwq6wEAAAD5MlGs5CVk8s6MFgC55/uyACgaoQVAroksisLntICRhBYAQItEFTAZoQUAkIgAg84q0ufDhBYAQMbsIgsTK1JgDRFaAAAAiQktAIAW2D0QslFvVivPM11CCwAAKKy8xpbQAgAACievgTVEaAEAALkzFFJ5D6pGhBYAheCobAAUidACAABKKcuD1wgtAArDrBZZc8RByK+87WIotAAAgFxqFE+tnp4FoQUAACWW15nYdqOo2d/L+n4LLQCALsp64w/oDqEFAACQmNACAOgAM1fQPXn6bNYQoQUAAJCY0AIAADKTx9moFIQWAABAYkILACCxkZ/P8lkt2KLezFVZZ7MihBYAAEByQgsAAMhUGWe2hBYAAEBiQgsAAOiYarV/eMaqjDNXjQgtAACgq3ohuIQWAECHDVYHHX0QJlG2+BJaAABAR5QtnlohtAAAOsQsFr2sUWT1SnwJLQCAhCaKK+EFvUNoAQAAuZDi6IR5eUNDaAEAAIWWl7gaSWgBAAAkJrQAAIApydsBLvIwwyW0AACAKUsVW3mLtnYJLQAAgMSEFgAAQGJCCwAASKIsu/2lILQAAICkBJfQAgAAEhJZWwgtAACAxIQWAECX1QayHgHQaUILAAAgMaEFAACQmNACAIAeNVgdzHr1XBVXAAADS0lEQVQIpSW0AACgA4py9L2yxFbe7ofQAgAASExoAQAAJCa0AACgg4qyCyFpCS0AAIDEhBYAACQ2NItlNqu7Wj0gxmB1sGMH0RBaAABQRgOVrEeQVN6OKjgZoQUAAAzLc9DkeWxjCS0AAIBIG3JCCwAAIDGhBQDQhLtn3Z31EKBjurFLXpF2+0tBaAEAAB09Al8RTXVZCC0AAJgih3FnLKEFAAAJ9PV9NeshjDNRAGY1ezU0czbV28/77JvQAgCAhMo2u5X3oGlHvfuU+n4KLQAAoK6J4qOdMEkVM0PXk+cI3CrrAQAAQFkUYTYrz3FSJma0AACAcUYGWZ7ibOxYGo2tNjDxZZq9f+3ed6EFQG6NfJEEoHWdCqRmr7eZy+Up4iLSjUdoAQAAozSKjW5+Lit1gHX7e8KEFgAAFFy7nw1rZc+BViMlq5mqsfdpKuOYyu8KLaCnjH3ytWsa0ClTeX7x3ES72gmuqa5vRdw9sFXtjF9oAQDkmOiiWZNFVpFjp4hjF1pAz7CxAs0r4kZN0XhOKpe8P2Y69b1TI6935Gegur088vh4ElpAz8njkzHkjcfJ5Cwj8qpa7Y9qtX94Hc3rutrMuKYabO3e9xShKLSAnlQbePrJ96wbsxwJlEve39Xvhrxu1NI5ma/3A5Wu3lwz97fZx0EWy66V25zK41loAYQNI2DqRr6BU+85xfMMnTTZocvzsP7lYQztamfsQgvgKUV+AQDSabSxOtHs90TPHyMDbCKeg4rN3+9pKWap2r2OVneX7ORjU2gBpdTshs1k1zH2+ryQAmN1+3nB81AxdG2XuCnsNpjndakMr7lCCyidet+V1e6TdRme6OvJ+wekKa5m16mVg4OdHEZdKTZ8UzxmOvX9Wh7P+VFvXRs6QAX1TfT47FS0NtrFN9Vrv9CCHKgNPL1LSlk37LNiWTY2NrbM2pHK2IjKIqpG3m4nbz+Pj5eVg4PD9znzgzT0qKH1olFwFUG3Xw9aua1Gl83b47FSq9Wav3Cl8nBE3N254QAAAORaX61W22WyC7UUWgAAAEzOroMAAACJCS0AAIDEhBYAAEBiQgsAACAxoQUAAJCY0AIAAEhMaAEAACQmtAAAABITWgAAAIn9PywG4oGcDJvcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(9, 9))\n",
    "plot.single(comps[4], ordering[4], fig)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAKACAYAAADn488NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3b+S21ieL/gf+1bFTHSkam5E34ro1hoqXSPLXEPtN/AkLW8eI8HHGGtVT0LIH3nrVDqSnKqO6JaxqowxtmcDa2QiE0QCJECCxB9+PhGqJIGDg0MQZCa+dc7BqiiKAAAAAGC+fjd2AwAAAAA4joAHAAAAYOYEPAAAAAAzJ+ABAAAAmDkBDwAAAMDMCXgAAAAAZk7AAwAAADBzAh4AAACAmRPwAAAAAMzcN30Kr36/KuJ/nqop7PVfEfH7sRtxwRz/afA+XC7v/bR4P6blEt6PS3iNzMNCz8Uffo34w9iNYDQf/nR/Dnz6U8SbX5+Wvfn1/ueQmuo8xX6W4M3LN/Hly5f49OnTP4qi+H5f+VVRFJ0rX71cFfHvR7WPY2wiIh27ERfM8Z8G78Pl8t5Pi/djWi7h/biE18g8LPRcvMkisrEbwaiyiFhnEUXWXmZVW9/0vLSrXFOdq5b1ZV1N68s6d23bVEe9vno9TWX37WtXHYcqborIsizW6/WHoij+vK98rx48AAAATN8QF5dcnl3hTtP6fc/77LPvvuvrdoVITXU0BTZN67vsu+s+uhyfrbI3+8tXCXgAAACAkzkk+DnXPvpud0yIdar6SyZZBgAAAJg5AQ8AAADAzBmiBQAAsEBN8/C0Df9oWQzMiIAHAABYhCHutLM055j7BJgGAQ8AADB7Xe/EMyeHBFNLev1APwIeAACAB1O8vXg1tGlrm2AHEPAAAADPTDHoONa+19QUknQJV4Zux75tAZoIeAAA4Iyqw26ODVGGnltmiqHOKdq0L7Q5VYginAFOyW3SAQCYjUPmWamWrW6/q64hLsT31VFvS7VN1eWH7K+pvq7b149NWx3nCCu6tOPQOrsuP7chXytwWfTgAQBgFpoCmbJ3R9cwpe3xkENz6nX1vViv91hpakfbPg4NBrpst+sYVdvV9bbch8wl07Suy/u/y5hhiiAHGJKABwCAWTrnXZMODRaGDlyODTJONe9LW4C2b5u2IKuPuYYkc203MF2GaAEAwAHmeIE+xeE/U2vPoY4dIgdwLD14AACYtCKLyMZuxMxNNWg4Zy+sczjVfE4AXejBAwAAF0gPk2E4hsBUCHgAAGDBBBCn5xgDUyDgAQCAC3PIpMgATJuABwAAAGDmBDwAAAA9nKPXk55VQF/uogUAANDRzfv7nwIYYGr04AEAAACYOQEPAAAAwMwJeAAAADowLAuYMgEPAAAAwMwJeAAAAABmTsADAAAAMHMCHgAAAICZE/AAAAAAzJyABwAAAGDmBDwAAAAAMyfgAQAAAJg5AQ8AAADAzAl4AAAAAGZOwAMAAAAwcwIeAAAAgJkT8AAAAADMnIAHAAAAYOYEPAAAAAAzJ+ABAAAAmDkBDwAAAMDMCXgAAAAAZk7AAwAAADBzAh4AAACAmRPwAAAAAMycgAcAAABg5gQ8AAAAADMn4AEAAACYOQEPAAAAwMwJeAAAAABmTsADAAAAMHMCHgAAAICZE/AAACxEkd3/m4optQUAlk7AAwAwczfv95c5V9hSD5nK503hkwAIAIYj4AFg8ors+QVseWHYdvFYXd9UX/Vxfftd67tcqDbV1aWdfV7HrvJt+29qW9vradpXUzvGcmhPlVO0v8v50KUdbefQMW1qOh/2nRO76tq1r6blbeXbnjedmwBAN6uiKLoXfrkq4t9P2Bp220REOnYjLpjjPw0t70ORRayy7ecR28vq5cv11W3Lx9WLinodTfvq87xNU7m2dh5aT1lXF23HsHqMqseqqY19j0XT8S9lD/+Yhiya34+mc6KpTMTz87LtYr7pM7nvwr/vPnbVUa/nGPva0LeNpSz6fT4O3U/9uNbf7+q6Y+ptkoXvAKYhi2Wei1ks83XRXRbOgckpisiyLNbr9YeiKP68r7iAZ04EDONy/EdXZPt/8Rx60dLVqesf2iEXXMdcpFXrGPo4ZeGPjinJwvsxJVmc5/3Y9f1w6u/HLJxzTEMWyzwXs1jm66K7LJwDk9Mz4PnmHG2akq7/BxyYp1OHL3MKdyLGG8Iyt+MEdOfzDQDTZA4eAAA62TdXEwAwnsUGPCbmAwAAAC5Fvzl4Vqu/R8Tn0zUHAAAAgIpXRVF8v69Qr4AHAAAAgOlZ7BAtAAAAgEsh4AEAAACYOQEPAAAAwMwJeAAAAABmTsADAAAAMHMCHgAAAICZE/AAAAAAzJyABwAAAGDmBDwAAAAAMyfgAQAAAJg5AQ8AAADAzAl4AAAAAGZOwAMAAAAwcwIeAAAAgJkT8AAAAADMnIAHAAAAYOa+6VP422+/Lf71X/81IiJ+vLs7SYNO7eerq9m0/eerq4h4OtZfIuIPI7anr6ZjXV3W5b2oH4M+++66XddzYm7Hf4l+vrqKf/7zn/Htt9+O3RRG4L2fFu/HtFzC+3EJr5F5WOq5ePc/7iJ+P3YrGMvVb1dx9z/u4ur/u4q7F3dPy17cxdVvV43bVNftKrdru+qyLvscSlsbTr3fvn788cf48uVLfPr06R9FUXy/r/yqKIrOlf95tSr+86jmjW+VRRTZ2K04TPbwj2GlSRKbPN9bLgvHf2xpksTHjx/j9evXYzeFEXjvp8X7MS2X8H5cwmtkHpZ6LuZFHpGO3QpGtYne50CSJ5En+bPHu8qdqo4mbeWbts2T/Flbmuqo76tpH/W2trVvn81mE1mWxXq9/lAUxZ/3le/VgweWqEu4AwAAc1JerMIpVUOMXWFK13XH1LErANqlbf9tdXXZV31d33Yd+lrMwQMAAADM2iGByFT3dWj9FxfwzHV4FgAAAECbiwt4AAAAOG+PB+D0zMEDAACwQF0mf318vDlLk4AT0oMHAABYhEMnaF2Sfa/zUo4DXCI9eAAAgNnbdcebuarfrrnvdsBlEfAAAACLcWy4McXbi9eDnqY2CnUAAQ8AAHAR9oU3TXPTVOesGTv4aZw7B+CBOXgAAGBkU7xYn2KbxnCK49ClTscf6EvAAwDAbBxy0TvWhfKuOWHKdcdOiFut6xi72jN20DDUa6zXt6veU73moV8LQJUhWgAAzEJ9eEp1uEzbEJr6xXT9eduQm2PmYWnbZ1vbmsKV6rCgrm2plqvX1afdbcOAmo5533Z1Xd61niHmy5lS4DKltgDzI+ABAGBW6kHEIUHPruW7QpJdgcIpb0/dte6utwk/9K5MuwKzMpTqE9r0NcUeRseYc9uB6RHwAAAwaXmSRxR71ndY1mt/O5Y1Tb7bp65D2zB1XYedNW3X5VguzaW9XuD0BDwAAHCAOV2gT7Wth4ZCU9T2Gsa+8xZwOUyyDAAAF2oJwcrUmVgZOBc9eAAAYMH2zT8kfDhO1+PnOAOnpgcPAABcGGEDwPIIeAAAAABmTsADAADQgZ5PwJQJeAAAAHoQ9ABTJOABAADo6PMPn8duAkAjd9ECAADoSS8eYGr04AEAAACYOQEPAAAAwMwJeAAAAABmTsADAAAAMHMCHgAAOjOxLABMk4AHAIBO8iSPIhPyAMAUCXgAAOikyJ4e50ku6AGACRHwAABcqD4BjTAHAKZNwAMAQC9lT55qj54IIRAAjOmbsRsAAMA05EkeSZ5sPW8j3AGAaRHwAADM3OcfPsfreN17uyKLWGX5IG0oA580GaQ6AKAnQ7QAAHjU1hOn3mMHAJgWAQ8AwIK0BTSHDKEqsqd/AMC0CXgAABagGuDUb2HeJdypBjnm0wGA+RHwAAAMrG+4cor9ls/bgp/yZ72Hjt46ADBPAh4AgIGdOyS5ed9vn/tCpyLrHkzVQyQAYBwCHgBgtuq9UaaoyCI2eX72fdb/1dd3qWNofUIjAKAfAQ8AMLpqALLJ88d/1XWb/L6nSPkzT/LH8KQtjKgHQPXeJtX97hpW1WW7qYcS9bBnqACn7yTMUz9OADBX34zdAADgstXDnV1ldoUTZeCzvfDhPw8/V0m5bcM+88ryPCJNKusfFpfLqtutkqd9r7KnNtRfS/k8TZLYZZPnkSbJVhCS5Lu3qVrCHDpFtn2sAYD9BDwAwFm1hTh5kldzl51lmww97Khp303L7oOd7vWXAU75s1y2yp6CnGqoFNEcdpTbRET8Je++/zGVxypP8k6hVddyAIAhWgDAgXYFNYcMw5l6OLHLoW2vDkXrWkd9m7kdt663Yp/L0DcAmAo9eACAg1V7oZSeAoc8IvYPSbpE+4aidV0+V+U5ssryiGgfglbtHQUA7CbgAQCO1hT0VNdFxNO8Mvm5WsVUPPbCyfesf1zw9FDIAwDdCHgAgKN0mST5cV37ahZq3xCr6hCzVTb8Xb4A4FIIeACAXjZ5HlnLcqjqE+40PQcAuhPwAAAwuGq4I7gBgNNzFy0AAE5GuAOH8dkB+tKDBwCAwbk4BYDz0oMHAAAAYOYEPAAADGbfxMoAwGkIeAAAAABmTsADAMBgzL3Dkjm/gSkT8AAAAADMnIAHAIBBmYcHAM7PbdIBABiUYSzTVH1fVtn2slUWAMycHjwAjG6T54//qs8ZlmPKXBVZxM3758FR9XmRPT3fVe4UqvuttqO6fl8bTtnGLm3atf+249pUps+6+nGrr+vavr6a6mpb1tTGm/e7tz/H+dZnH33Lwyl0OQfPfZ4e+tk4tJ3neH2roii6F16t/h4Rn0/XHAAAAAAqXhVF8f2+Qr0CHgAAAACmxxAtAAAAgJkT8AAAAADMnIAHAAAAYOYEPAAAAAAzJ+ABAAAAmDkBDwAAAMDMCXgAAAAAZk7AAwAAADBzAh4AAACAmRPwAAAAAMycgAcAAABg5gQ8AAAAADMn4AEAAACYOQEPAAAAwMwJeAAAAABmTsADAAAAMHPf9Cm8+v2qiP95//jNr8M25MOfhq3zQ7yJiIg38WFr+c9XV/Hj3V2/ug5s24c/3f/ssm3TPurbf4mIP/RvRu+2DLXtrtdU11Z3vXxZrqy7rb5dddfb1fX9Peb4M4wPf4qI/4qI34/dEkZxqvf+1/vfF/GnD/vL7CvXZT/H1nMqv77p16b/ioj/p7LNr7XHEffPy+VNx7B+TKrb7FpWeX519XPc3f34fN2ufexS3Vf1ePSpp+k1NNR3dfVz3L2421rWum213oYy//Iv/3d8++233dt4gLu7H+Pq6udBy5bvXVm2bbv7cn4TM77775z/FVdX/3hctu98bzrPq8932fp+e9aOHxs/O+Xjtm2b6ok4z/cI09F8fvienZo3byK+fPkSnz59+kdRFN/vK78qiqJz5X9erYoP2f3jIjuwhSNLkyQ2eT52Mw6SPfy7JKvs/mff822Vdd+ma9ksLu/4T80qi4hNRKQjN4RxeO+nZYz3IysistXT44jt5+XjHZI8iTzfNJet118v07aPrIgkSSNP8u39lM/rbW2z6/Xt8er/ehWvX7/uVDYiIk/ySPKkc/ku2/Wts1q+PFatdeeb8JuYKUiSNPL8L5Ek7w+u4/58vq+ra9l6G/rsp7pddVm9no8fP/b6HmFZfM9OU1FEZFkW6/X6Q1EUf95XvlcPnoj5BjvM06HnW5/tnNMAM1ENO56FL92CkDzJI5KWsrvq37GPerjzuJ++7du3/z0ODW361lcPZqpl6q97VxjUZb+P5dqLw9m9evUuIg4PQroGNNWyeb7ptV25bX27vnUA89I74Jm7PMn9kQAADGZXWBERB4U1Q6n3Kqou29vulu3ry/fV02d9n7bBJTk0mBHowGW5uIAHAGDpmgKSY0OTQ7Y/NvwBALq7uIDHcBwAYGk+//A5IgQmAHDJ3CYdAAAAYOYEPAAAAAAzJ+ABAAAAmDkBDwAAAMDMCXgAAAAAZk7AAwAAABcsSdJ49erd2M3gSBd3m3QAAACgWZKkz5bl+WaEltCXgAcAAABoDHfOVa8Q6XgCHgAAAKDVqYIfhiXgAQAAABbrkICq3qOorGPKPY0EPAAAAMCoqiFMGaK0BTNNIUuSpJHnm8efu7bv2565EPAAAAAAk7EvXGlbXy4/dThTD6OmMjG1gAcAAADgALvCpK5B01BhkIAHAAAAoIOuoU2fXkT1socGPr87aCsAAAAATuKQYWYCHgAAAICJOHQOIQEPAAAAwMwJeAAAAABmTsADAAAAMHMCHgAAAICZE/AAAAAAzJyABwAAAGDmBDwAAAAAMyfgAQAAZi9J0rGbADAqAQ8AALAIQh7gkgl4AADorMjGbgEA0ETAAwBAJ8Id5ipJUr17gMUT8AAAcBCBD1NVDXQEO8ClEPAAANCLYIe5EO4Al0TAAwBAb0IepuqUoU6e5CerG+BYAh4AgAs1ZEhTZEIfABiTgAcAgL2EN3Dv8w+fI+K+N48ePcCUfDN2AwAAGF8Z4Kyy58sAgOnTgwcA4ILVh1YJdbhkXefv0XMHmCIBDwDAzN28F8zAUKZy5y0hEtCXgAcAYGBTD1v2tW/q7YdTq4c8SZJOJvgBaGMOHgBgtopse86YKWqa2+Zc+6zr0wYhD3RT9rRJ8mTUdgAIeACA0ZVBTT1UaFpWLq9uW19Wr3fffnfV1basb2BzrqBnVzBzirl2BEHMVdkjJ883O9cPJU9yIRBwUgIeAGBUhwQETdu01XPIcKS2ZW13mCqDqGogdWjwUd1P31BoKWHLHHpmzUmSpK0hxiVqGn6V55tBAp2m3jzm0gHORcADAJxVnyBmaoHFvrZ3ae+u3kr79tPWlg7FJ2GM4WrQxZDhTv0xwLkIeACAg7RdrLuIP0zXMGtqodch9NBhLFOaKLkaAhm6BQzBXbQAgJMosuOHKy1V3+NRHsslHceuvZSW9JrryrBhSqED48iT/PFf+RygLwEPAHCUfRfh1ZBnyRfrHK7p3FjyuVK95XZTyFOuHzr4qe+DaRLuAIcyRAsAGMSSL8g5nUsKdvZpCl3qd3qqPu8zebLeQtPglurAKQl4AIBeimw+k/oybcKd7uphTpfApk/44y5b53VML52mbQVGQISABwCAExLaDKdv75su5afSo6d+q3KBU7O2YEjPICBCwAMAwAm4UxZ9Nc0RVA19yuf1stXlS1JkEWly/7hrj588yYU8cMEEPAAAnITeO9PQdNGfJ3kkkT4+jvx87enTa2jXvESH1jknhwzlqvfm0bsHLoe7aAEAAI263s3rXAHLUoOcc6jehh1YJj14AADgQjVO2HtgiCJ8mZ6m97ct5Hn18dVpGwOcnIAHAAAWrm2YVmlfOFOdBLnpbl4AjM8QLQAABmPeneWq3pp96GBnLkOHnN/AlOnBAwDAtGWriE1Esnr/uGiJd006taFDlCFCnmqbxrwDlImIgSXQgwcAgNmZ0tCgskdLvWfLoW0sJ8Othx9t66ag3tYu7ewzP0zXfR9Sfqh2TM0SXgPQjx48AAAM6tKGsdR7ntRDnrbeRuV21Z9t5ao/m/bbtm21XUUWscq215fvVZrE0bqEObteZ1M9u+YN6rrv6jHuUs/Ofbbv5qya3ksAAQ8Ao9p1Idj0x+u5/qgt23WqfY3xx7kLAuaqyCKyiHifjNuOLdlq62me5JWwJNlankT6rGz98SbPY5XEM13CkF3Diqrtimj/zm2rpxo+1b9DqmHNru+Xptdb3Wbftl1CoV2qAVmSJ61BV7ms+p5ERKyS7eN28z5iXeucVT+up/yubXoPq8vKfbe919W2Vcv88GrcYXKc11Dvddd6Nnm+9d14Coe+pnO07VxWRVF0L7xa/T0iPp+uOQAAAABUvCqK4vt9hXoFPAAAAABMj0mWAQAAAGZOwAMAAAAwcwIeAAAAgJkT8AAAAADMnIAHAAAAYOYEPAAAAAAzJ+ABAAAAmDkBDwAAAMDMCXgAAAAAZk7AAwAAADBzAh4AAACAmRPwAAAAAMycgAcAAABg5gQ8AAAAADMn4AEAAACYOQEPAAAAwMx906fwavW/iogf7p/86cMJmjOc6xf3P29/215+9dtV3L24O3+D9vn1zf5j+l8R8fsJtOOc+2krN2Q7u9Z1juPPbr++iYgvEfGHxtVvYtrfS0vx4U/3P9/82r6+bd2u+rY8fC6r9ZTv/Ic/3a8v3+8P8eaxzCHnwG/X9z9vb9+0lrm+fqr39vZ+39X9dlGtY9/+5qH9s8gYnt6Pq6uf4+7ux7i6+nmrxN3dj1vPy/X15dX11XXX1x+2ztv6fsqy1edt69radX39IX755epZ+1++vIuvXyO++67t9cN53N6+iT/+8UPc3V09Lmv6vFXd3f249fkpfx/88svVs3L1z0z9d0fZhqurn+Ply7vHOpo+Yy9f7r7mKT9r5T4+f/6X+Pbbb3duw3JU3/uS79npefHiTXz58iU+ffr0j6Iovt9XflUURefKVy9XRfz7Ue0bXZInkSf52M04zCYi0rEbcWZZEZGtprGPSzz+E1NkEau4iSLWYzeFE1tl9z+L7GlZ9vBvrlZR+31b+d6pvs6+dRax2np+PlnM+x1ZllevfojXr19Hnm8iSfb/smoqt2/bcn213NbjJI/IivsySR5JnmxtGxFb9VfLlH+bPT6vteXmJo937yLevt370uCk1uskPn78GK9fv+68TdP5XNa1q1y1bClNi8bP3679pOn974YkSbfqK/dffj7L7xEuk+/ZaUqSIrIsi/V6/aEoij/vK9+rBw+c3anDnXPtg8HcCHcuwqGBx5RVg5iIGCQbqdf5bB81qyzuQ+0OZXfWc9YgiXN5FvjUQpoyuIm8fZvHMh3qbpMneSSRbpW72dd4mLD6+V8PdtrKVcuW4U21TNPnrypNi9hsVlvL2/YNLMPFBTyz7b0DAEe6D86E2ktVDUbK51X1sOZxm8ryLurbNv1tVWQRq+xp/VZvnkr5+ra71sEl69Izr2mbLoFOWffHj713AUzMxQU8AABLt6+HTNvyemBTfdwU0uwKaLo879tGAKCdgAcA4EI0BSdde820hTqH7hcAGJbbpAMAzNznHz4/Pt7VQwcAWC4BDwAAg1niJOkAMAcCHgCAhdBLBwAul4AHAAAAYOYEPAAAAAAzJ+ABAACAC7ZeJ5Hnr8ZuBkdym3QAAAZXZBGrbOxWANDHep00Lr+5yY+ut6yjuo9j62WbgAcAgJMQ8gAsQxnK1AOZanDTtY62ZcKe4wl4AAA4GSEPwHLsC2nmYqnBkoAHAIBRFdn2c4EQAEPaFUL1DaimHAgJeAAAOKkywKkGN/VQp2kbQQ/A5egTtDTN54OABwCACSuyiHwzdisAmJKxg50u+6/29EnTIjab1eO2p+oFJOABAOCs9vXe6VsOAM6la7hULZck6dbzeh1DBT4CHgAAzkZoAwDbhgp8BDwAAAAAE1EGPknSb7vfDd4SAABooPcOAJyOgAcAgGf6zJMzdHAjCAKA/gQ8AAALVAYv9QCmfCxEAYBlEfAAACxANbhpC2+6hjvCHwCYHwEPAMBC7ApmhDYAsGwCHgAAAICZE/AAAMzczfvDt22bowcAmBcBDwAAjXP3NE3SvGt7GFOaFmM3AWBUAh4AgAsllGFp1utk7CYAjEbAAwDAXnrxMHVJko7dBIBRCXgAAOhNmMOcrNeJIVzA4gl4AADopAx16j9hDvTwAZbum7EbAAAAMKRyLp4830SEcAe4DHrwAAAAiyTYAS6JgAcAgM4My2Lq6nfSGjLkyZN8sLoAhibgAQAAAJg5AQ8AAEBHn3/4HBERmzyPTZ6P2xiACgEPANBPtor4y3rsVgCM4ub92C0AaCbgAQB6KTIXOMAy1efvaaPnDjBFAh4AAIAHu0Kec04ybkJnoC8BDwAAQEU95EnTonPvHoCxCHgAgEEUsbqfnwdgAdK0iDQtImL/rdb1tgGmQMADAAymyOI+5Mkawh7hD4yiDCnYVgY4bccnSdJIknTv8dvk+VmHbgG0+WbsBgAAl6HIIlbZKiJzsQmMpxxqVfbKyfNNRNwHPk09dfb13ilVe/EUWUSaJI/Py0mZq8sAhibgAQBO4jHQiXgMdYosQj8eYGrawp0+6r143GkLODdDtACAQRXZ9oXO47AtYIuhU+fXNFFyORTr3KoBkDl8gCEIeACAw2Wr+8mVW5RBj/kp4LkkSWO9TmK9TrbCnvLx0HdtEihNzybPt+bw0esHOIYhWgDAYbLVU4Bj4BV01jQcqD6Zb7VMmhax2Rz3GSvrru7j2DoZVhnulD9XyWhNAWZKwAMAHESvnBHUJ6l+nONolNZwoH3Dgarr1+vk2fOqm5s8IoYJgTiPTZ6bbBk4CQEPAMBMbE1czeLsCn7ahmtV7wi1XiePd4QqbTarSJI08nzz+LOq7NEzpXBoyWFV2TsnT/KIfHfZ+ztxPV9ezteT5MnWYwABDwDAGKq9cWp3G6urDoEre06tstM1jfmqh0T1YKh+a/Cp2Dc/UD2ImmIw1ccxPSDvg956PXoFAQIeAICzK2J1H9k09cbZE/Y81pFVNhmoXSxX05w/Y2sLderzEdXvclWdn2joiainpuyhUw+E7p/nz8qXw78MA4PLJOABADiXypw5Tf8Hf9VU9qE8DG2su2rt2u++4Klpguol69PTpz5JM3B5BDwAAGe064KtiNVWmGMYFmMre8iUkzm3qQ6ZantcnS8oYnuYWNceRpcU7gyl7AVUup/bJxmjKcCJCXgAAE4pez5/TlfuVMa5rNdJ3NzkrUOe9k18XE7gXL+9e7muqV5hzenlSd74PfJ4K/bMBM2wJAIeAAC4UNU7a+2bz6Y+tKraK6esq14349oXElfn8vnh1avTNgY4OQEPAMDQqnfIeqA3zuHug4Us3r5dj92Uxdjk+eOcT12CmKYyS5/gGGBuBDwAACdQ3tp8FcKdodQDhX3zwpzTviFMQymH1phDBYC6343dAACAJRPuLF85dClNi0jT4jGIOqaHyybPtybH3eT51t2R6s/31dVUdgpDqIa481OfYxHxfNLheluObQ/AWC6uB0+RzeeOFOUfhHNp76kc8p45dgDTVUR7L4dVjHPbZjhU22276yHPzU1+t3dsAAAgAElEQVS+NelwxHaI0NQjpzo/yj5devY0Tbg7VI+gTZ4/1lEPR9Ik2fla6yHPrrbsCl7a1pX7L+utH9dq+6rl6vIk7/p2TIKQiinJk9yE3mfQK+B582vEf2b3j4e+cJ5y8HJo26YUMhzTlkO3bTpu1bqqf2C01V2WaWtDl4njmrapLpvyuQdwaYpYtX9HZ0Xj3DZb204kINJr57TKi/AyPMnzzVbIcuqhW+WuVnEf7tRDgeqFTLWtTcFO1SEX5PVtdtXRFu5EtF98VV9bU5kuvV7a2rjKojEwaQpZhgormuppa181cGq7G9XSbPLcnbXoZdf3X1WRRaTPF5+0LV31DZ/q3wdTGjK7KorufwitVqu/R8Tn0zUHAAAAgIpXRVF8v69Qr4AHAAAAgOkxyTIAAADAzAl4AAAAAGZOwAMAAAAwcwIeAAAAgJkT8AAAAADMnIAHAAAAYOYEPAAAAAAzJ+ABAAAAmDkBDwAAAMDMCXgAAAAAZk7AAwAAADBzAh4AAACAmRPwAAAAAMycgAcAAABg5gQ8AAAAADMn4AEAAACYuW/6FP7222+L//2//zsiIl7cnqQ9g/kQbyIi4k182Fr+89VV/Hh3d5J9/nZ9HS9uhzswP19dRUQ8tvdLRPxhsNpPr+l41I//vmP223XEL79cxcuXLyMiOh/f8thFxN73+7frbufzvuP/2/V1RHRv49QNfT7X647of6x+u474+jXi7u4qXr58/r7+8stVw1ZMWfnZ/uWXXxqXV339+jXuTvD9fXf3Y1xd/fx4TtXPo7u7H7eeX1393FjPy5d3O8/Bly/v4vb2zePz6+sPrWXn4OvXiO++qy+9fvi5jO/B6bt+fPT169f47vkbsiC38fVrxN/+9mZ/0RHUP8/Vz/q56ri+/hC3t2+26inr6PJ9s29/TfW2lanvt16+7/df12NRHoM+y7vsq/o6rq9/i9vb/zeurv7xuL78/r//3XUbEddbv9eqvx+qvwuqv0/KbavlIq7j9vbFs/aU29V/77x8+fLZfqvfE/W67n/3lW2O+Pr1j/G31d+e7Y9lun7xcL7/6eFz8OubOOUV5/1n5/n5PEfnfC1v3kR8+fIlPn369I+iKL7fV35VFEXnyv+8WhUvkiQiIjZ5fmgbR7XKIorsNHXnm00kaXqayiMie/g3F0McjzRJTn6u5ZuIpEMzs5jX8Z+yfLOJiOh9fqRJEh8/fozXr1+follMyM3NzbNl7969i8+fP4/QmmW5uckbl6dpEUnty7Bedr1OHpe/exfx9u3gzaO3++/Td+/exdtFvyFpvHsX8dNP3f9unbPNZhVpOtxr3WxWe8sMub+umto1Vju67nezyePt27dH/S2S5w9/B+35A/Tm5ibSNHl8Xi+f55utZTc3N7Fer5/tpzzO63XyuKys7+Ymf3ztf/3r2/jp9U8HvCIWISvCFU83m02+9dk8paKIyLIs1uv1h6Io/ryvfK8ePBHzDXaYp3Ocb13CHeDc0igvXBlWGdLUNV1otJUFTutUIUe13i6hz6mNEeY0OXc79gU7Xcs+X5dGRPKsXPn6kiR9ts16nVSWvercLmCazMEDwOQIFqCr9OEfcOnqvzvL4KYp2Gny0ye9d2DuLi7gOdXwrIj+w00AaFftZr5eryPXgxQuVpoW8dNPz4du0l29l8pUes9wOnmSj90EWKRzDc86RO8hWgBwLtWQB5iPNE1is8nHbgZ7lCHPFIZrcRp9hoMB83dxPXgAADgtF5XTorfOkJzbwHQJeAAA4IIJgLpb0hxxeZIbxgULY4gWJ2NOIto4NwAATiNP8kjypHF5RESSJ1vBzuNjN6+E2RPwAAAALEAZ1jQGOC3PgeUwRAsAAGBGuoY0whw6y0y2vgR68AAAAExQfVhV6/AqGMJf1hHvs6ewJzM/19wIeAAAACasy9ArGES1J8/W4wHDnnq9AqXBCHgAAADOoNoLp2kiZJisliFc9V5l92Ubgpq2IWCGhg1KwAMAwGDSNIkkWY/dDJisXT1vqutu4ub0jYEjNZ7PUwxtdrVpQb2IBDwAAAwmSdKxmwCzYZgVVDQNCasHM6cIYNqGpbWWn24IJOABAAA4M+EO7DDnIV3ZarQQSMADAABwpLbA5vC5dvSGg8naFzTVQ5695YcJhAQ8AAAAR+g6r04f63USoZcPzFff3kYD3LXsdwdtBQAAgKFWwPGawqADhqMJeAAAACZGcAQTdY55gA7ch4AHAAAAYOYEPAAADG69Xsd6vR67GQBwMQQ8AAAAPRg+BUyRu2gBAHBi5e2eN6O2AoYk5AGmRg8eAABOJk2TsZsAg8mTPD7/8HnsZgA0EvAAAHAySZJGmhaRpsXYTQGARTNECwCAk0qScojWzdbysndPuf7mZns9ANCdHjwAAJzder2OJEkr4c994OPOWwBwGAEPAACT8BT2pDvLQZN6YAhwaQQ8AADAIrizFXDJBDwAAEzKep2M3QRmaFe4o3cPcAlMsgwAACxKGebk+eb+Zxn+5OO0B+Ac9OABAAAWSa8d4JLowQMAwORs300rGasZzFCSpI89dpI8ichWo7YH4Fz04AEAABbJpMvAJRHwAABwVtu9c2BYQh3gUgl4AIBe1usk8vzV2M1gpoQ7AHAaAh4AAM7EhLdMm0mZgTkT8AAAcBbrdTJ2E2AvIQ8wVwIeAIDFcGEKx8iTPPIkF/IAsyTgAQBYEHPcwGkIfYCpE/AAAFyMtPYTaJMk6VOok63cnQuYvG/GbgAAsAzrdRI3N/nWPCs3N/lo7blEu+e4qYc7aURsTtoemKIkSSPPd5/7W2FO3lYKYFr04AEABlMPGKrP1+vEJLtnkz4bqpWmRUQ8vQ/eCy5Z2TunOuyqdQhWtjpTqwCOowcPAHBSgoRz+Cki3j4+K495NeRJkrX3gotXDXHKXjpJ7O/RAzAHAh4A4GBlYGAo1rh++ukm3r4duxUwfU3z6ORJHol5qYAFEPAAAMzcq1fvIiIZuRX0tdWbRA+S0ZlEGZg7AQ8AcJD6/DpAN0mSRp7kj3P3JnlyP89LVmyVKQl/Tk+4AyyBgAcAYCaqdyqrD4urT6rMdNXDhMe5YFom+S2X5/mmtUy5rikM6nLXKADmT8ADADAjekstVzX4SfLk2fOmcKc6UXDEdhhUfV5/LPABWB4BDwDACJqCmrbJqpvKNvXiOU558e/Cfwq6DBmqlnkMevLkfkG2egyJymXVx3NU7YlUD7sEVgACHgCAs+na+6Z/uY+HNOehjvuhXTc3ZWiUhpBnerrOEdMU+pR3icpj+rcGL+cnqs5HFNkqIuJ+zqJ8x3YnbRnA9Al4AAAmZKwhWOV+73+Woc/NKG3h8jwbfvYQ6rSVqw9fA0DAAwBwcl2HU5073OkyMfNTDx9hz5z1uUtUU4iy1aOmZZtdPYK25v8p682KZ9sdcjcrd8ACuCfgAQA4oe2eMfPizlzL9xiwVHrM5EfU1SRP8ojKPECPslVEZfLovvMOAbBNwAMAwMFO0cNnk+eRJslg9TGAbNXYi6frkKldwYzQBmAYAh4AgBMoe+wIKwbWEjTQz32vmu0JmPdqmBenabs+gY1wB2A4vxu7AQDAsmzyPDZ5PnYzRlUfqXLO4zHWBfNZhnOVAUO2ap2El+7ahlQBME8CHgDgJPoEPet1Mss5avo6ddDTN9wZuj19Qp5Fh4BlADXxEKp6G3Wmp8jGbgEwN4ZoATAZ5QXfJQxnGWPYzlD73OR5ZBGRff4cEfvfr/r7uivI6Xq3qTkqj0P5PjQFHLuWN9VVLn/7+XO8f/16a/0h7/Xu7dKIaL9LUpM8ySPJk8cAoT43S57kWxexc/rsl5MT750geEIhT5FFrLKxWzFvRRaRRcQ6fb78XMe2/Mx4L4G6VVF0H8O8Wq3+HhGfT9ccAAAAACpeFUXx/b5CvQIeAAAAAKbHHDwAAAAAMyfgAQAAAJg5AQ8AAADAzAl4AAAAAGZOwAMAAAAwcwIeAAAAgJkT8AAAAADMnIAHAAAAYOYEPAAAAAAzJ+ABAAAAmDkBDwAAAMDMCXgAAAAAZk7AAwAAADBzAh4AAACAmRPwAAAAAMycgAcAAABg5r7pU/jf/u3fij9+/eP9k+v7Hy9ub4duU2+/XV8/X1g269mq26aFwzhh1RERX79+jf/jb39rXd94HKpO3L5e+6ufNvva1fp+HtmOHr5+/Rrffffd8RVxxPt5G1+/RrS/Ddfx7A1/3Nf431VExO2eN71c3fB2fY2v8d11+3fg8bp+YRFxff+d+Lfvtj9b5fvb+Lv3Yf3ez+J13P72VOb6xXaF1XVTdv3rde04VL+bql+CDa+n03GKrc/T1/ga38Vyf0fdxouI+BIRf4iIiOvrD3F7+2a09lxff7hv1+2brcdP63+L29sXj2XrbS23qarW9bSw4Vx4/Jz1/CzsOa/ajue+Y30dvz3U23I+79S+TZdj1lV5bBuP8d72RWx9fm8jvv7xa3z3XcPvo+rvuF3vz235/VD7/ry+3X6fWuu7fmjS7fb50LZtvY56uYfHX2//+PB79jp++eWX9vazCHd3PzYsffqe7arr9/F1/PbwXd5t+SEO/d3Qtw3X8dvW8y7bVn8v9PHmTcSXL1/i06dP/yiK4vt95VdFUXSu/Mcffyz+4z/+IyIikjTt3bhTyzebDqXSiOhS7gAnrDoiIk/TyHatf3j9SZp2PBYnduLjcW7v3r2Lt2/fjt2MZSi/PnqfH2m8exfhbVia+omQRqTPT4538S7ebn46T5PYY1P5Tqz9PZBuGj7b22XS9w+1/OV5zeW6udv8ZRPbX3bVX4r15QeqfE7exbt4G28Pr2vi0kgiInv4txybzerxcZoWret2Kbcry1fr2WxWkabF1rr6Ppu2O0S5r6a2N+23rY372tL1uDTtf9+xaiu/VSbdxLu/vou3bxt+H6WbiE35We/52X7cds+yp9Y830dbHY+bNNRV2eZd+tf4nHzu1WyWI883scTv2VPYRP7we+n0iiIiy7JYr9cfiqL4877yvXrwTF2X0ClNkri5OUNjRlB9/ZMJeQD2mt7/MGCf8j1ruMDZNAUZLbUsJMzZL6383DQsP9CmOQxlPrqEKm2hSZe6di0bKthp2lfXtgxVpimMOaR9++pJ0yI2kUfjd19EJUA54LPdGuS0tmZvHfXgbJj9AlO2qIBnSPWwaKlhyb5Q7P51pw9ll3scgIlx4cpiuZiim3rQcEwwM1SIcwpl245p47lCqqXLX+XxOl6P3QzgCBcX8GzyPCLPT1T57tVlmNIUkqTv04fu3NPSFABVlwl8gLPYpBHvxm4E3QkxmoalNQ1JO1r5f9/Tv56gcqboEoMHDpMkvovh0lxcwDMFTaHJw/+7eHy+iiI2m/xcTepsV4+fZ72eljQBD3B+evEwY5cz/Iy5EhQBHOZc8+8cQsAzUUWsIn94XIbveUvZPnb1IjqmPujFNTtdmRsA5klAC7NRhn2vXv0wbkOAowl4BlQPO1bZ/c9Dh14N26tye2JHwQwAnNchw7R23XEMGMEmZjVkuOlOaeVk0uUQrvLnx4/jtBEYjoDnhIqsfDROmPI8xDlzO/yPOwB4ZlfQY2gXMLQy2EmSNPJ8Y3geLJiApyc9X86tPN670qJ0q7eTiZ8BmKoywClDHoEOzM15rwXqd0/rYr1O4uYmf3xenWy5DHlMwAzLJODpQbgzhC6BTaX0+7InlGMPwLIId4Cuqr1uqkOuqsHPep08e9wU5Ah3YLkEPB0Jdw7TftzSyDcRZdDTNvnz0zA3ALhc1bl4qsGQuXlg2cowJs+f/8/RalBTDXfgEEmSxsePH+Pz52zspnCEXgHP7W+3kb6//yIxcpNj3f9OSmvLxgjStiegbmqHYV8ATIFePzBv5fCpeu+bfdrCHMEOQ2vr4dUUMjI9vQKeN79G/Gd2opbA4NKHOQ52fxndDwPTQwuA+SkDn7+O2wygB5Mcc6n2DQ8UIh3PEC0WK0nve5rlmz09g7JztQhgWerDhgwXAuim6UK37N1T9sq5OXObYBdzN82DgIeTmcq8RfvaYZ4fgON0GTYkALogm2n8/oepahpuVd71ypArOI1dAVXbneXqPYp2zYk1FQIeFidJywmcAZgC88YAPGcuHdjWp5fQrpClb2+jtvJz7LUk4GGRZvhZBLgIwh5gycphVuXjUrmsl00aIfiBRlMKX6Y0MbWABwA4iXqYI9wBlqwMdJp646RpMakLUuA4XT7PQ/VI6kPAAwAMSpADXJIuw6uEO8Au9e+IQwMfAQ8AMAjBDsBwzMsDl+sp8OkX9Pxu+KYAAAAALM+Ue+QJeAAAOA23TAeAsxHwAAAwPOEOC2b4FDBF5uABAHpJ30fEp4h4PXJDAAB4JOABAOB0qj150mFuAwtj+umnm0iS93rxAJNjiBYAAEBHr169G7sJAI0EPAAAjKPs3bOJvneChVHotQNMmYAHAIDxbFoeAwC9mIMHAIDzewxz0ufLzdUDAL3pwQMAwAh23EbdLdY5QJI4b4DLJuABAAAWQcgDXDIBDwAAsHjCH2DpBDzQQZL6gwAAYI6SJBXuABdBwAMAwPRsUnPxAEAP7qIFAMB0bdKItBi7FcyMHjvAJdKDBwAAAGDmBDwAAJyXoVeckN47wKUS8AAAcD7CHQA4CQEPAADnsRm7AQCwXAIeAADORO8dADgVAQ8AwFIY/gRHSZK0cQ4f8/oAcyDgAQBYEsOgAOAiCXgAABalpafBJn3q4VN9DOxU9t7RiweYum/GbgAAAGck2IG9hDnAHAl4AACWxjAt2Ksa4uS5Dw0wfwIeAIC5++tPET+9rSzQ+4AnSZIKMCr0zgGWyhw8AAAwMqHDebQd57a7ZwHMiYAHAABGIFAAYEiGaAEAwBlVg522x6cYUmWoFsCy6cEDHf3w00+RpGkkqf/bBsAEuTvW4pTDhvb19GlabygSwOUR8AAAwMR0DWHKctWffQKfKZpTWwGmxBAtOFCSppFvdHMGYOYee/48/E5zbT0LXUKcLkFJWWZqQ7eq7XI7c4BuBDxwhHK4lqAHgPl7uIjeRETq99rUDN2rZcz5eNoCm7n3PAIYm4AHDmAeHgAWoe36vuzVk26eyvjVtziPvWSSPCLfU+aEYVDfEKdeXq8egHvm4AEAoNkmjftk57zpTvr+rLtbvDzJ95apz+FTn8unz3Cvrsu71ruvvImjAe7pwQOMpt4Tqmmom7mOAOA06sFP37t1NQ2vOiTkAWAYAh5gssoA6GmuozFbAzB963USNzd5jy0GuOiuDucaSJH5zj+1Lr169hHaAEyLgAeYDHMbAUtVZBGr7HT1//TTTbx9u95ZJk2L2GxWp2vEnn1HxOH7z1YRm4hkdT92y5wr/ZWBTpIno7YDgNMR8MAAmoKJfLN5HF60L7gwBAk4h+q8Jpu/jNcOTuAv9+FOnuQRvXvxDGSTtvbiOcUdoCKmFfQ03RlqzDtVtRmi585Q8iSPJE9GD5/KdkzNVNsFTJeAB06kPrzoElVfezXEOuaY/PQp4u0RbYK6c/VqKHswPMqe9nlo2FK2PX1/X0efiWmHmsT21D1TlqDIupcpj2V9m33H+KbL+5ltnysnCfk229/vj+f9w/l+Uys+Zq+iISVJuhWcJLWhb12CnurF/BQv7LuEMF3bXT1WbY9LQxyHLu3aVabLa9/keaTJ7n3sU29DkUWkO6oc6jwpX9+rj6+OrgsYl4AHJmDOEwl3DWuGCLrKC9K2C1M9Es6jPP7l8a6/H/UL00P+53XbxWf5P8irF7tDvO9pWjwLXNrOs4MDjWz7IrZaz8Fhy8MFe9c6hg5j9gUX5f66hBXVtlXLLzE8anut1fW7tq1qOz5Nn6Fyv+n7p8dt53p128fP/I4eOlv7rgzHqgcb1XYlaUQRq4j06Xuib+Azek+eh891vqNI0+TD1Z4r1Qv0etjRJ/Cp94Z53H8tMCjVz522gKJeX1v95fK2HkKHBhFN+2s7HnmSP4Qizfsqtytfa9Nr63Kc60FI9TbzTcexrafSrtdx/6D9NdRt8nzr+7Z8fdX9Vtc92xcX59xB8hAhKPutiqLYX6osvFr9PSI+n645AAAAAFS8Kori+32FegU8AAAAAEzP78ZuAAAAAADHEfAAAAAAzJyABwAAAGDmBDwAAAAAMyfgAQAAAJg5AQ8AAADAzAl4AAAAAGZOwAMAAAAwcwIeAAAAgJkT8AAAAADMnIAHAAAAYOYEPAAAAAAzJ+ABAAAAmDkBDwAAAMDMCXgAAAAAZk7AAwAAADBz3/Qp/O233xb//d//Z0REXF39fJIGndrLly/jl19+OVHdd/HLL1cnqTsi4p///Gd8++23B29/6vb13d/Ll3ePj9vKvbx7Gb9c/fJYtk/7X969jNt4sfdcLfexz7HHf5e7F3f7C03U1W/9z6lD3s+IiLu7HyPiS0T8Ia7jt4iIuI0Xj49fxO1j2d/iulOdL+L2WdnqsvJxWXd1eXVf+54fss+m9bvKtW23rw3VMm372bV+l11tLB/f/uk2rn+9X1Z9HNcPx/D2/vnX+BrfxXd798l5VN+P23gREfH4WSxVP59NquvLx2Vdc1d93X2Oz228iOvrD3F7+2Zrm3odZbny83H7x69x/bfdn4/b66ePVb8X03Gj2+bvhIP3W6/669eI73wHMK7r24fPW/VcLM/v6+bnt5Vt47r2PLbL1rdv+uhU67mu77uyfctH8lldZbk/fv0a3/mMXZRn55fv2cl58+JFfPnyJT59+vSPoii+31d+VRRF58pfvHhR3N3d/2GRJOnhrRzRzc1NrNfrE9Wdx3qdnKTuiIiPHz/G69evT1b/0IY4Hnm+Ofm51nUfpzz+eZKfpN6p2vzl/mff82OT57GKmyhiHau4/+4qYhV5bCIiIomn97Fctk8SaeSxefxZXVZfX623vq99z5v22bZs1/pd5dq2a2pHvZ5d7a6vr7djl2q5dbKOTZ7vrGur3s1DW9L75+/iXbyNt3v3yXlU3480koiI2ER+/749vGdpJPfLWpTblTaRP1s2Z+Vrr76m6vGoHp/y8b7XXy9Tbt/l85Funj5WTxU8vV/tO+34e3hfPUdINxHx7l3E27cn2wd0sUkfPm+bt08L04itX4nlR2ZTe/6w7HF1bfmWdLvs3vL1bcuHlXK1X6v3VVTKvvvru3jrM3ax0gjfsxNUJElkWRbr9fpDURR/3le+Vw+eiPkGO6VThTucxhzPt6awJsmTs7djqW7i/jNcxGrklsxXEod9rsrt0kiiy/8aqO+nDHdYlqZAYl9I0bZ+SeFOxGlezzF1ds1pjtIlMIILt6k+2PO5bMqIDt5f+bweSAGL0TvggaUZIkTa1wPn0nrocLhDw5dzagvXuvbq6WSTRpJG5P4AnZ003Y7/tnrxbNJn6y9VGdRMOtTqG9acJUGCGXiWqER7QtP2Edu0PO3ze3FTC3M6BErAvJlkeUCnHJ4FlyBP8q1/Q5hSYHKutpT7OfX+qj16qhepUzrmnFjPC3rhzowIa6BZU8DSJXTZdCw3him3DehFwANMVp/AZ9//BR8ydDikrvo2TXWMEYwMsc8iVscPmdukWxeUZRAww1GaF0Vg09+oPXaENnAZdgU2Tb2LgMUwRAs4u7K325BD13YFDNUeLbuGEA06xKhh//XHY9i1/6YQaojj0eU114MccyzBCJqGZAmF4Lm/jt2AIwl1YLH04AEuSteAJXkYeNRl+2NDm2O2b2vnuRy9bxePXKBJTzLtMwkAsyXg4WSWPCeRSZOPM+QcO4caOxg5xrF3wTrHvjrXP8+3gANMIrygmzFCntowTaDBIT1v9NaBi2KIFhxo7ICC4ewLMcYKgupDpHa1I83SKLLxh4CVut5KnXm6D2vykVtBRDyFIu8i4qe3+8v1XQcAzIYePMDsNfUMmHMPnYj+QU2RHbf90MyhQ5XeO/QmdAKA3gQ8wOwNGSaMHYzU9ZkL6JC6Z8GFHnRT/6z47ADQ0SbmP384hmhBZ59/+Byfk89jN4MZmE1wAsxeeav6v/61R9DddLcsAIj2aZuO/et2U6mjug9/NQ9LwAMwEkEQQAebNMwUCzCupm/h8tu561+0TXUIe4Yl4AFgMEOEViZonrFNGpF692ZBLx4AjrSp/ZyTpQZLAh4Azm5XEGSC5vlKhTsAwATtCqH6BlRTDoQEPABMWpK6ITf0Vk6wbKJlAGaiT9DSNJ/PuZTD0sp9T+k3rYAHAGAp9v2lK/ABYAHGHha2aXlcldbKpA2PhybgAS6OyY0BYrywp+v8P+YJAmDG6r/BdoVCQ/1GFvAAACyGABsA5qYp8Dnkf3H8boC2AABAP4aLAUCjQ/uvCngAAAAAZk7AAwDAeTT12tGTBwAGIeABAGAc7voFAIMR8AAAMBIBDgAMRcADAMB56ZkDAIMT8AAAMB3CHwA4iIAHAAAAYOYEPAAAwOzp/AVcOgEPdJAn+dhNAABgDyEPcMkEPAAADO+QK21X5xxjM3YDAMb1zdgNAABgwTZpuPLmbFpOtTI7TJ2KwIIJeAAAODE9czivaqDz2DFMuAMsnIAHAABYJB3IgEtiDh4AAGAxHkMdwQ5wYQQ8AACMKn1//w8GJ+QBLoiABwAAWA6hDnChBDwAAIyqyO7/wZwkaRpJagJxYDoEPAAAnIaeFMzMZm9ec19gvV6fvC0AfQl4AAA4kcN6N6RpMXA7oLv9IU/EJs9P3o5DPz/A5XKbdAAAgIjHXmeb2I5XHkMfvdKACdODBwCASUkSPRcY3yYV7ADzIuABAIAF6zLk6BKVAc6+41MPeX769NNJ2wVwKAEPAABwMdpCnWfLNpV/FTfvn29r0mVgCgQ8AMBJGGYDTFolwBmml1Mq6AFGJeBhkfIkH7sJABepiFVEdv/PdzEwNcfOqdN2h7fyrlrVu2slaRpJKugGzkfAw+K4oCQN7k4AACAASURBVAB4Lo1kuMqyVcRf9v9f6iIbbpdclvPcgnp8ZdiwSbfzhsHnzDFB8JOGIVe9ZKtni9I0iYiIp06LaeP63dKO5QDaCXhYvKbARwgEXJoinl+UnLvuJEkbL444j7aeB4MbOJ1Y4ki/6hwwbYerS5lD9zlUnYtyRPBT/w5M38eznjtJ2i3ASd97Y4DDCXi4CHmSb/2rLgOWJ6n839MkDvtjOc0O/yP70H2eQxrJ0b15iuxpktFynp3qBU6RVXrvVAKdTZ7r1TOiqc2JtF4nYzdhXA0T+DZN/tv0vOnuT50Cm023u0bRrst3WL1M+f1YD4KqIVCS3q/3HQkc45s+he9e3J2qHTC4PuFNU9kkTwZrCzA/RRaRj92IiiTSyI2z4EGSpJHn/c6HTZ5HfnOiBtFP/a3bNCxPH5536OnT9DgiIt11itSCpZ1lGUyaJpEkadzctH8YkzSNfOMNgf+/vbsPkpo8/AD+jSi+4KlYsdx18PhhC9KK1rGFEzuzeaitghaqcIzSKlSn2ilTrUWh9e2ydYr4Qqe1jlVsLTgdp0UUR8u0VjFZLYxIGat1sMWKRdADrC9wiMhbfn/sJpdkk2ySTTbJ7vczA7ebzT55Ni93u999Xii8UAEP0P9BmB9+KQnOoKUR55lXEMTWPUStI+nwpFb5cWxfQIYaIJIyvkG2bi/oc93UGitFV5Bg57DWpclathLICDRZQw/KLXnkvL+YJBmXqlerG6/HrTNDqS7ru99NTahwyrJ+3oIpHRKgAZKsQUd1V65+tZpZCQihQ1W18j0hm7eJqHWFDnhalSZrDLUiiCMk4b4noqxqRHATlA4p8rbcnssWQw2gSIASbFwcI5gTsgxNU6ErgIjQiidLdAXQ1NYZUDk0r0NrDXR8gpuq5QFb8zSaa+sh1R7gmA+7rJvboEfxf1wW5evDizC6yQqjm2w587WGPkTUehjwhMDWS+nz6krF1jZElAShCKhK+R22Ma5O3KFHmPF6vOrQqDF/nC19ZFnkvhVJVji7XHl1wSqHIc3TJspvWKBaH3Bbjle3rjDPyYha4wyVF9p+2Je5rK86A6+mUX6x5etBLQ/UrBhj9VT/Lqgn7Cei/GPAE4FXmMDgJ37WAZFrrUNErcFoWRK1hUmY52VtHB6nJFrZGAMwu7UrcX5waIZuQmny6r5mDIZca5wdVdMyHfUwoMm+QBlRJTipFSPbWtxYVjbu25Y7xxoKWhnUKKMJOUNQHRKguK8rhO4xq2HInUxEucVZtBLgDBzSmq3JbeYoIiLKhnpm6fIt12eWrCCzZ+mQXD8guD03zGwvWZvBqdFqvX6jm5IOCZqslQdEljXIsgi075Ys8R89uZ7dL4TeuCnWqWGMKbuFkMvjt1gHa7b882s146aqoVGQqd7D5g5urZmYXRARMeCJm3MK7nqCFeeU3nHUrZFThDNUIqKkZHka8qDqmQrX7/V7BTTGYzYh5kr2KjMIWRYtPc6KX0DjFuAY50aYc6Szc3Hgda1hjTGOhx+/cyoJLZ4FJsIIc3yXG1OoexVSGQvHddr2yuM1xwFy3na7XwuDHE9e16lbSCsL4TPAMxHlFbtoZZBbC6Aky4+7LHZVI6KsCdOVybluFgYb9gt06gm7GjWEDrtyVVgGVTY+iNUKvuoJxry6SOmQIAkdqppgaKMKz24z7LqVYQHH9XEbF6dmmTzmoXl3uQrHHKS9KPtOz05E+ceAp4GaoUVLkBmtmuF1ElHGqAJYDGDJrJqrhg086m0NZA2A4mpZlIVQyRDmNXm2HFJiqkxO2f52hpg5K4xZmzaFWt+tVZGwhD7Gw0GCmFrrygKQKqM6hQ2VzGoKGYAGdYl9TBfKKB6bWEQdyL5YLKKnR4PzQKiaBgEw5CGqW3ZHdW/JgCdMANHqrVHc9hUDHCJqFlEDmXoGeAbin4nLa1tZCYlakiPIUTUN0ABJ8X5KI2clUzUNWk//VOUGkVDLHiP802qsV7N7ltF6pPIzKx1MVMD+fj+DAZTRHcuYQluHZAZvVE0IOXPTjReLxQAjqfX/vnELXcvXvj3gMc6NmTMX11dBIkpd0wU8QVqYhC0viXWzIo91JiKqVyPCjygtiZKqk4Bc9THObVmaas0YlUtKOdTITQsmJctzcmVbVRhgneI7Y2GP8WHeuP7jvvbCBCPGus7wKQuyMN24s4tmfGOZCdtt41yYZfmfiPKpaQZZdhtA2O1fMw/yGPegzFkR5Y1xbt5ME1Gi4uwyVe92otQlrvq7dZ0yx2QI9H1wMqz1Mv8+N1HIoCvVf490xXvQ5TjeoxgDqXrNeuU3G5bb385iUQ703FakC1EeYNiYHrwyG5UZVljWtQ5D4xUZWAc89hoUuV5VA3pDCnXeedXLOliv8/eN23Os+8n5ut3++W27nnq71SfK41lVrnP17xsOZE5xKhaLAdfkidcIoQKeM3v7b2dp+u0s1CEJzjeGm4YH61/vtj/qCTzc3qDWU1aQZXFh0EPUHGQBDF8Sc5kt8EZDQK75Ov0G8HR7btL7zQhAjH86pEp3J/d/1nXzItYvm4x9gf4PbaIE76DMstwIb6wz7OiKPcix1jXIYK+1ZuWK/TAJGXCZ4hslURm3J5paE8yZoQ7kSmu4yjGonIui5KxQpUuZIyhwhgbOQMNWTv9acPugZKxvC0dK1TMleV0v1g9p1hmWhJBRLBYhhFx+bVXhS2W9yratr0OI8r7x/ACo9O8318eU8gxuRjm1Zn0y6ur3uPFanQGS8dqKxaJZhihZjqUieRyPfHObXct6LJ2Cf5inrKv6/VMSrsudnOeA82+a1zli/cLAUROf5SFbQ0ecGS7atS1CzUYnSqIh14+k68G/jZEk6V0A4UbxIyIiIiIiIiKiqDp1XR9Sa6VQAQ8REREREREREWVP04zBQ0RERERERETUqhjwEBERERERERHlHAMeIiIiIiIiIqKcY8BDRERERERERJRzDHiIiIiIiIiIiHKOAQ8RERERERERUc4x4CEiIiIiIiIiyjkGPEREREREREREOceAh4iIiIiIiIgo5xjwEBERERERERHlHAMeIiIiIiIiIqKcY8BDRERERERERJRzDHiIiIiIiIiIiHKOAQ8RERERERERUc4x4CEiIiIiIiIiyrlDw6wsHXusjqFDk6pLQ7T39aG3rS3takSzcydwzDFp16J1NXr/921o3Lbq0TayYZsauQHYiZ3YOvIYjLTsng0jYbufpEZuK6s2VA65137Y4HNK1LPvjGPvV05fe/l3fJjtuB3TvvY+tPW2ud73Ogecz3Gu29feZy5v622z7Se/+rqVG1Wc56/b8XCWb9x3nhNudbA+17lvNsDy+kf2ARsc973WTdFI9NfLec147RdjnbDLAWDD0J0YudX+N8p5PJznsfX5XpxltKH6CX2wV8w419t622y33co23pdZrxPj+nDWfffu3TjqqKP8K+zQVnnP19fXZ1vW19dn/vRb16vMWutQc9N7ddt7EevfHuPctV5vzuvAet/r94Nxu73P/rfDYL1WDNbfE0YdANg++xjlGden9TrV9V4cs9WowEgAxk+roH9ErM/fEGB5FkWpo7G/3J5nLa/2fujra0dbW2/I7QfX19du3i7//tuAbduOxYEDB9DR0VFZx/670rjfX0YfOnZ1lM+ltl5bmd7b7UNHxy7LkpF45513bOUa6/XXzb7cuSwsZ13rLdNtv1j/1hh629ps17SxbrlO1esDQEdHB9577z3897///Z+u60Nq1UXSdT1wxaVRo3Tcf3/g9SlmixcDs2alXYvWleT+L4lkyk1bQY21OFUAi7EYszAr1nIpH3jss8U4HkItX5t+ROVXQa31AEBArrtuWaFCA+D/+oXj16Qqqpe5PddZppi5GOqSWb71KSklFJSCf6V9yHA/gBqqK2zdltd2hQooJe/HrK+5pJSgqiqECPf3slCo1KFUsi0rlUrmT791vcqstQ41r4JScP17FOR3oaGklMyy/BjXSJB1g27H7Xo01nvzTcXjra5xjTfp+9UWViopKBQUlI+xwOLFwKZNMnp6tMoaamW98jmiaRpkWTafr2kaeno0s5xSSam5TeM5VsWibJZjXQ+AbXtudYjCWdd6y3TbL7IsVy1320fW1+lWj56eHiiKgmKxuE7X9S/VqkuoFjxERKHUCq5iDoCIqPGCfKCpCidcQhwVWlOFO0G57b8o+zQIpVBwiWLiJSBDhWb7AOn1wbT8GqIHTkRZEuaaDBPW1BPKFpSCGd74lWcse3OmV0kMdppVOdwBvI+xqKwH3/DGCMgB42+NtTy1Uk75L5ARaGSNM/z3CvPrCfmLRRl1ZlM1cQweoriVRPU/v8eatfUOEVEIzRzuRAljwpq5JP16yBDQISW7EaIWUvc1q4q6AiKi4IwgSHGEO/2PlX/6n9T9gZN1WcFcbg2SrI/73fdaVk1U1rPW116v8s9sf3ZjCx6iONQKaRjiEBG1HCO0UqE1JOQhIiJKSqGgQNNkWwijaXJC2ypUtfTxC2mCdA+zhzdez7H/sbZu0y18qkdS3X3ZgoeIssOvZRNDMiLKKSPoyWIrJaHU/7tVhvAcn4eI6tCIZJi95SkmxaIcSznOIMatRU15PcWjRY1wCWbcWw/1l+FWh3DKLY3iazEXtSwGPESUniAtn9iNjYio5bDFE1GMVPgEObzYKGv8zsn+rmBO9mWizsDFv1uZe2see92s3cqqu5FVPz+ucIhdtIgoXQxviKgFWFvxGLNr5YEMAQnlGVeDz7va/NwGryXKrsp7LVW1jnVrWV75ubixtSKKmzPkadT2+ruquW3TWKZUnlOwLFc8yi041g2OLXiIiIiIMqLR3bh0xb5t45+125UOiYMnO7FbC2VRzWn5gg92S0TxcpuxzC2Q8utSFgQDHqIo2OqEiIgSYm3tY/zzW7feUMhWhsIwhyj3VOHdz5H9H4liGS8o7kGX3YW/XtlFiygqhjxERBSBNZDx6q7VyJY8RqAjKQ3bJBGFZW01JtAf1AhHczIGOEQtjQEPERERUUqCBjnWKdfDPjcoa3ctIsoIM7BR0T9mjtvjRETsokUU3H+XcEYnIiJKXRzdsogoJZEDGb7/pHyLaxp18seAh4iIiCJhyNB4fvtc5gdAovxgyxtqccWizNAnAeyiRURERHXL2/TfRESp8ZvtSqgMf4goMrbgISIiorqwJU/6lmB42lVoGoVCIe0qUCtxznjFcIcyIEzrGrbCyRYGPERERBSZW7hjHSPG+ZOS0dm52HZfKPyQSNRwRjgTNKRhmEOUO1kPtNhFi4iIiBLBcKdxZm3aVL6hCvPNp3NMHg2q62xcRBQj19Y4KsdIJqpDsSijp0dLuxq5wBY8REREFBu/MIezPzWO25Tnxr7XIfFYEEUVqXWOKE9tzhY71CJqdfGyPuZ1O21u9Wpk/aJuiwEPERERxYKBQfqE0D0f0yEBimTe1iE1qlpEzS1QcMNwh1pPo4KROMsPUlaWgignBjxERERETUTVNM/H3Fr2EFEM2DqHyPdLhlrSDk0avf2ktseAh4iIiKhJsFVOPDiTVj4VFB43oqjiCBxkuXFBp1c3sKDLkhK1BVBcdWTAQ0RERESNoaZdAWp2ikc4x/CHKLwoU6U7nxMluAgz/k3Qx4MsrydkieN11yozCAY8REREFIqAjCUYnnY1KIdEid1YKEGqf08ptswiCi/trlNxi/p6vAKhWoNE1xpwut71nRjwEBEREVFD6Er19O15wGAgRxytxMyWO3G2Hpu5hGPuUEurN4RIWtjWPI1Qz/hEYRzakK0QERERERE1iC3YKSW0EYY81GQ0rRP/93/Vy4tFGT09mnk7qrgClrgCnFrrxxkIGeMTJT0tPAMeIiIiImooGQIaB+Sh2Am4NdVhCyyi/PALPcJ2dWpF7KJFRERERC2p1gd/BgM5poKDehORqVXCIQY8RERERJQ5MoT5L2l+QQ5DnrzxO19qn0tes3AREVllNQRiwENERERERE3DfUr0YEEhh9UhcpfVQIPsGPAQBcFpXYmIiBpGQE67CpRnKlAOdIz3byHex7FbFxHlGAMeIiIiIiJqUvySjohaBwMeIiIiIsoUHVLaVUicezei9MppJWHH2THGYeJ4TESUdQx4iIiIiKjhhFJpWaHCHE4ZqmV5ypwf5guFQrwf8OPqCsQuRS78zyGvcXaMsKygFBicEVEuHZp2BYiIiIiouRlj6uiWZboCCLkIVWiAUl4mC/s6WVcoFFAqldKuBkVgDXBKiuMYqmDPLiLKJQY8RERERJQoWa58Wtbsy1WtvEBX4tsWQxeqydHqqSAK7iEPTyMiyhl20SIiIiIiSkEjugGxq5Gdazc71fGTiCinGPAQERERUWLkjPd1SXXgXNU/gHE+FiasMddtktCi1nGKL8hyP185wDIR5QEDHiIiIiJKRmU0W1XTzO5YrSBUGBAggAkb7BSUQi6CnXyFJtkOKomIAI7BQ0REREQJEJDzkDFkihl4GFmCcxwYtTxeDOAyMLBlnUDbUvrHnbHedqtTPWMa+ZVdT1lG6GWWHevJJgColgBKOH42hhA6VFWK9Nzx43sxcODQUM+RZfv9vXu3YvXq9kjbJ6J0MOBJSe/48Rg6cGC4Jzl+627duxftq1fHV6kWEnr/yzrQY1+0dddWtC/kHz0iIiI3OiTvz8O9vcDQ4B8+ewD09Nj/EO/atQsLFy6MXsEsc87iFGDAX7cgxWjJE2fIkqoEBj4eP348Bg4cCNl8n+02j5t9WaOCEHNw8gjChjtJlUFEdo8++iiOP/74UM9RFAWKopyJ/l9G2wC4XqAMeFISOtxJqIxWFcv+P5p/9IiiGt87HgOHhrsO5co0y4a9W/didTtDbqJcChHueDn66KN9H89X9x8XzpmejG5XleVVMz+5hR+WMnxDHo/gxNiHfq14arbwSSCUsXVBq7PsgTG8J2QQQkRBhQ13PHza6wGOwUNERA0XNtxJqgwioixIfKYrS9CjFApQEgi/GjpbV4DuWMWijGJRtt2345g6RNR8GPAQERERZUhvL6Dr4f4pzgW9vWm/DIqbM9SwzMBltrTxC1kqz1eFOfZ1YvUzB3puMCF0lyCnvNxyD2mNqUNElDQGPEREREQZEkPvpZgKyQe3rlg1p9TOWvetuAYIDlKOCltXr7D89l3Ux2qt61uuJUyqZ8waIqJmwDF4iIiIiIgaqDyGTHmmproEfnoM2wrDMT5OqHDHMcV7VdjjnALeOSA1UNU1ywh+6pmVKg2qpkHrqb0ekZuwM6k5Bw8HgPfffx9Tp06Nr1KUOAY8GbNnzx709vaivb0dRxxxRNrVySWj2bHg3KyZ1zunN9xg1T0caJcoCWEHvS53dlBsy7ZuBdo5sSClpN6pvBvOfI/iHrwIIUNVtUBFFYvFqhnGqsvTPctrdGumqoAmLLfnBizP2sKnWJTR06PZHj9w4AAGDBhQ9bzt27fjxBNPDFHJbPniF7+I008/HV1dXejq6sLpp5+OQw5hR45mF8fg3zENCEw+du3a5TppwL///W+MGjUqdHm8sjNi+/btmDx5Ms477zzccMMNOPfcczF58mRs37497apFV0qnmaxQBIRi2XaAevT29uKDDz4AAKxatQpPPvkk9u/fn1QVqSKOmcg40C5R/eK4jlqoRxDVIEcY1+Tmm2/Gxx9/bFu2efNmXHLJJXFVKzVuAUqxWLS0MinvL03TYB8fJn71Dq5cbxhUKAQIdxL+gs45+LLVN7/5TXzyySe2ZevXr8e5556bbKUStm7dOvzwhz+EJEm4++67MW7cOAgh8OMf/zjtqhG1vIkTJ+K9996zLVNVFRdffHGk8tiCJyPmzp2Lq6++Guecc465bOXKlbj++uuxZMmSFGuWX6oI1orn+uuvx0svvYT9+/ejs7MTbW1tOO644/D73/8ef/zjH5OvKFGLKxaLkKT+JvNHHnkkxowZg/POOy/FWhFRI40ePRoTJkzANddcg8mTJ2PBggV4/vnnMX/+/Ni2kXRLlThbEZVbmmiVe0qgVjrVyq2DisUiUGn9ysbN/mbOnIlJkyZh+fLlOOaYY7By5Upcd911ePjhh9OuWl0GDBiAM844AyeccIL57+WXX8Zrr72WdtWIWt6CBQtw/vnn45FHHsGwYcOwZMkSLFq0CH/+858jlRcu4NnQBgi5fDtgs1EKZsuWLbZwBwC++tWvYsGCBSnVqEkYrXecrXjk/tkUVq9ejVWrVuHAgQMYM2YM1q9fDwAQggP1ETVCV1eX7f6ePXvw9NNP409/+hPuueeelGpFRI00Y8YMTJkyBeeccw6uuOIK3HTTTVDVbMcR5sxVdQY75a5CAoBszvYky8IcL6ZQKFRa9xjrFwGgKvARojzWjKaV91sWdl/d3bEswnRbi2ratGkYPHgwJk2ahKlTp2L58uX461//iiFDhiS63aRNnToVkiTh5JNPxtixY3H11Vdj2LBhaVeLiACcffbZeOCBBzBt2jR0dXVh8+bNWLlyZeThWsIFPO3rAM66mQhd112XHzx4sME1aR62blo+Dj/8cADlbzeGWvoYWFsUEFFy3Jq+T5kyhSErUYWzlZvVLbfc0uDaJOOll17C3LlzMWHCBNx2220oFovYt28f5s6d21JjElrHiOm/XQ5//AnIlRFSqwcUliv3ywGJqHxZaw1Lyo8h/nGMIoY7QshQFHtdyqGXGmvII4QO61vwe++9FwDw2c9+FrfeeituvPFGPPLIIwCA73//+7Ftt9HGjx+Pl19+GW+++Sb27NmDTz75BOPGjcPJJ5+cdtWIWt71118PSZIwbNgwPPTQQ7jiiivMv+133HFH6PLCd9FSjDcYGfhaoIls3rwZY8eONYMeSZKg6zp27tyZcs0SEnB8HlVREx8secuWLea+37Fjh3m7afd9Bt177725fuNE8duzZw/27t2bdjWIMsHaym3evHm44447PL8Yyqv58+fjgQcewPDhwwEAsizjd7/7HQqFAtasWZNu5RLkNRaM+7rFyqDA1ud4hz/uU4YLGMOku4UoUUXt/uZVB1kWKBSsLZTKr7P8msrLvVoyhVEur/+N5qBBgyr1ErYvGfL+pd+cOXPM2++++y6efPJJTJs2DZs3b8b//ve/FGtGjWKEl274HjxdF1xwgXn7Bz/4Qd3lcQyeoEoCKCSXNGzYsCGxslNV50DLAjJUofW3xjGOgZD7w8Y6j0vT7vscWbZsGf+4tLDu7m7bm+c9e/Zgy5Yt+OlPf5pirVrXl7/8ZdvxsH7x8OKLL6ZVrZZmbeW2YMECfP3rX0+xNskwWklYfec730F3d3cKtckOt5me+tV+j2UESNauW0bwY9z3f74RomhwzpyXpnJ3tKJ5P9oYRdWGDRuGCRMmVC2//fbb6y47TWvXrsWaNWvwwgsvYP369Rg8eDAmTpyIcePGpV01ahAjvHTKe3jZDFatWoWzzjoLY8eO9TxOYTDgCaPkCBlqLQ/hoYce8nzssssui1xu7imS/e2Lsa8Vl2URcd+nb8uWLZ7fLDD4aX533XWX7f4RRxyBT3/60ynVhtauXZt2FchHs74Znzx5Mp544gkAwJVXXolFixaZy5999tk0q9YQfi15wrTyicLosmXwGlMoSCud4K2Cyu/dyt3KZJfHZEtXsupz3t46qby+EUbV49Zbb8WOHTtw4YUXAigPlXDVVVe5Tp2eBv/Az9v999+Prq4uzJs3D6eeemrT/h4hbzNnzqxadvDgQTz66KMp1IasOjs78fjjj+PGG2/E/v37ccYZZ6CrqwtdXV0YPXp06PI4TXoUJWH/F/Q5Pj7++GPz3/z587Fnzx7zfl31TJGupLr5wBLZ9xTKgAEDcPTRR2PQoEFV/6j5rVixAp2dnejs7MS6devMcCdKv2Oq34oVK8wuqm+88QZmzJiBSy65hLOtpKi7uxvTp09Hd3c3Xn31VfP29OnT065abPr6+szbr7/+unm72bqixcka/AQJgdy7bLnTNM0xsHPt8ovFYuBtCKGbA0q7PQaU62sMNu3cjv1+7boFtWLFCjz44IP47W9/i76+PkyaNAmf+9zncN9994UqJ8lQzro/vPahU7FYxKWXXooxY8bgk08+wQMPPID7778fu3fvTqqalDEvvfQSpk+fjunTp+PVV1/Fb37zG3R1deEf//hH2lVred/61rfwy1/+EqtXr8Zzzz2HcePG4e6778app54aqTy24ElChGDlqquuMm//4Q9/wJVXXplMXaytjBLudpYXie17Cmzo0KFsLdXCrF307rnnHlx00UUAgL/85S+YO3dumlVrSbfeeiteeOEFAOVv/H72s5/hhBNOwOzZs1uiJUUWOVu5NSNriwKv23mQ9FTszuAgbJAghG6GMOWf4Z4ffDv2cmVZoKenxwxnrN2rnJwhkSyLxFsxGY466igsX74cM2fOxB133IH58+dj6tSpDdl2UMb+6OnRqvZV/8Dadt3d3easdJdffjlGjhyJE088Ed/+9rfx2GOPNaTelK7Zs2dj4cKF+OCDDzBhwgTcfPPNWL16NQ49lHFA2p588kmsWbMGGzduBACccsopuO222zB27NhI5fGIximmFjOJvplJuVVP1uXtjWRema27Kt3lZ8yYgXfeeQdvv/02RowYgU996lNpVY2o5R155JEAgO3bt+Ojjz5K/AMr1fbRRx+lXYXEGRMeAMCHH35o3t6xY0cs5fM8rq3WODbGNOzV61QHMGFaC2WJMQbZwYMHIUkSbr/9dnP8nTjHIDP2l6aproEMUJnhC5K5bpCuWV77feDAgTj88MOxe/durF27Fg8//DCA8hcs1BqOOOIInHXWWQCAESNGxDKYL8XjxhtvxJAhQ9Dd3Y1x48bhtNNOq6tbKAOejDCmR9N1Hf/5z39s31pnuptCEwRGud33TeTtopiRaAAAB/BJREFUt9/GBRdcgFNOOQX//Oc/ce211+Lyyy9Pu1rUIMZ1Z70GdV3HG2+8kXbVWpIkSfj1r3+NF1980ewCtG/fPnZbTdGdd95p/p1yfhHx4IMPplQrf3KAAYCtrrnmmrq3mcUQJ2t1cgYAitLjutz/uQLRBlxOrsVQXGbNmpVo+dbuZ8ZPY5kz6JFlAWi1y3QGQG6teHRdx2uvvYZnn30WX/va18zl7KKVb2HGZDK69+q6jo0bN5q3JUnC0qVLk60o+XrllVfw0UcfYe3atXj66afx85//HPv27cNnPvMZLFy4MHR5LRfwqALlabdjGBjZSVcASfFZwRqGyPY+s8OHD4eu6xgzZow5Vdqrr76amUHdTE0Q6DjlZt83saeeegrr1q2DJEn4+OOPMXHiRAY8LWTevHmu1+BPfvKTlGvWmnbs2AFd11EoFHDppZcCAHp7e3HTTTelXLOcUSRAiWf8mG984xvo6OhAV1cXpkyZgnfeeQcAcO2118ZSfhbcfPPNGDlyJC688EK0t7cDgGugRfGrHrAY8AtiikUZsmy/3yySPA/L+01UzV4WtLVTrf1snTHN6Z577kFPTw8GDRqEO++8E0B5jDXr9MytIExrqLypFfacdtppuOGGGzB48ODGVaqJqZoGYf1FWKcDBw5g//792Lt3Lz788EO89dZb2Lt3b6SyWi7gAWAPKWoFFg0ao+bxxx/H0qVLbRfdaaedhosvvhjf+973ktuw8/UrOqBq3o83odT2PZkOHDhgG8B19+7dWL9+PQDg85//fFrVogbhNZgtxx57bNXsdSeddBJOOumk4IUolQ9DMQUcre5Xv/oVnnnmGQDAzp07sXbtWhw4cACTJk3CjBkzUq5dPLZt24ZnnnkGy5cvxwsvvAAhBKZNm4aOjo60q+bKrTuTEDJU63uoHIo6S1PQsv0UlAKgCmia/3pJSuo8rB77yF+wQa2913GGcGPGjMHSpUuxceNGPPjgg1ixYgVGjRrlOiV8GEmeL3FKup552A+LFi3C3LlzccIJJ2D+/Pk4/vjj064SVZx55pkYNGgQxo4di3HjxuGyyy7DSSedFKn1DtCMAY9XIGEENaoAgszeWKu8mO3bt68qUR08eHC45C6OQZMVKdz+aQKx7Huqyxe+8AXbIKKnnnqqeT+r3Q8oPrwGE6ai/KW8Wj3wqcEaw/z973+vGtjP+AY70BgUilR9v9WDHmMf1BF8Ga1Kb7nlFvP+/v37Y6ti2g477DBMnDgREydOxHPPPYdrr70WGzduxC9+8YuG16XWWDRW1lCn/MG9xyxDjvHb3VYgICPtqT/qOQ+tgYsma4DlQ3+Y2cW8xuRx0mQNsiMM82pVcPvtt2PlypUYMWIELr74YqxcuRKLFi0KtB238CJooGGsV08ZYddNm1croUa9But5aD0VRowYgWXLlmHOnDkYPXo0Ojs7w/1tp8SsWbPGdbDrp556CnPmzAldXnMEPEFCGLNLVrJViUqSJLz77rsYMmSIuWzbtm3hpwf12Rc1u5ClwBhst1a9gq4XRWz7niLbtWuX6/IBAwbgRz/6Eb773e9i9OjRns8XlabkapDO6pS68tgc/ddX6tdgJQARzi4JlQ/lqqqZAQmA/h4MXs/LGuH4WcOZZ55pzrYSlg4JUCy/q51hT6up9foDhl+SJGH79u048cQTIUT5QPb29uLgwYNx1DITSqUSHnvsMbz11ls466yzsHTpUpx88smBnhsmkAn6XGN5rbKtoQ7QH6L6zRIVhteHeACxdg9Ik6L0mK3HrS1d0lDPeQj0Byy6Anjk6Z6Mlx2k9Y4ma6HKfuKJJ9DZ2YkpU6bgK1/5SizDEHjVM0iQEaZbn7GuNTixbqPW9pzPD8MvoApbRhzbsq4T1r/+9S9cd911GDVqFF5//XUcc8wxocugZMQ9k5kU5g20JEnvAtgUaw1alK7rZ1rvP//885g9ezamTp2Kjo4ObNmyBY8//jjuvfdenH322Z7lSJK0LvHKNiHr/o+67wHu/3oZx2HTJvdfK7quY9OmTZg3b545bbMbHof84TWYLdbjIYSIHPDweMTD+R5h1apVmD17Ni666CK0t7eb18h9991nzoriJsvHw/kaDznkEHzpS19CZ2cnAPusln4DgGb5NVI+WM/FqOdhZd1MnYvOa2zLli1YtmwZVq5ciVdeeQV33303CoUCjjvuON9ysva6KDjrOXDBBRfgrrvuwimnnBK6HJ4D8XJem93d3VXjfOm6jr/97W/o7e31K8r1W6RQAQ/FqmrHv//++1ixYgV6e3vR0dGB888/P8hAWC3+9Whktv0fcd8D3P/1CvQL6KmnnsK5557rtwqPQ/7wGswW83js27cPhx12WNRyeDziUfM9wqRJk4KMoZDl42F7jV5BPwDzw7aHLL9GygfzXKzjPASydy56vsd6++23sWzZMjz22GMolWqOzZC110XBxfVBn+dAvBL9+8eAJz284NLF/Z8NPA6ti8c+W3g8sqUVjkcrvEbKh2Y9F5v1dVFwPAeyKdHjckhMhRMRERERERERUUoY8KRnW0bKaFXc/9nA49C6eOyzhccjW1rheLTCa6R8aNZzsVlfFwXHcyCbEj0u7KJFRERERERERJRzbMFDRERERERERJRzDHiIiIiIiIiIiHKOAQ8RERERERERUc4x4CEiIiIiIiIiyjkGPEREREREREREOff/d1OqyZTwPBkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x648 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(16, 9))\n",
    "plot.stacked(comps, ordering[7], fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q files?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Log-likelihood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
